A model for redistributing heat over the surface of irradiated spider companions
G. Voisin, M. R. Kennedy, R. P. Breton, C. J. Clark, D. Mata-Sánchez
MMNRAS , 1–11 (2020) Preprint 15 June 2020 Compiled using MNRAS L A TEX style file v3.0
A model for redistributing heat over the surface ofirradiated spider companions
Guillaume Voisin , (cid:63) , M. R. Kennedy , R. P. Breton , C. J. Clark , D. Mata-S´anchez Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester M19 9PL, UK LUTH, Observatoire de Paris, PSL Research University, 5 Place Jules Janssen, 92195 Meudon, France
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Spider pulsars are binary systems containing an energetic millisecond pulsar that in-tensely irradiates a closely orbiting low-mass companion. Modelling their companion’soptical light curves is essential to the study of the orbital properties of the binary,including the determination of the pulsar mass, characterising the pulsar wind andthe star itself. We aim to generalise the traditional direct heating model of irradi-ation, whereby energy deposited by the pulsar wind into the stellar envelope is lo-cally re-emitted, by introducing heat redistribution via diffusion and convection withinthe outer stellar envelope. We approximate the irradiated stellar envelope as a two-dimensional shell. This allows us to propose an effective equation of energy conser-vation that can be solved at a reduced computational cost. We then implement thismodel in the
Icarus software and use evidence sampling to determine the most likelyconvection and diffusion laws for the light curve of the redback companion of PSRJ2215+5135. Redistribution effects concentrate near the terminator line of pulsar ir-radiation, and can create apparent hot and cold spots. Among the models tested forPSR J2215+5135, we find that all models with heat redistribution are more likelythan symmetric direct heating. The best-fitting redistribution model involves diffu-sion together with a uniformly rotating envelope. However, we caution that all modelsstill present serious systematic effects, and that prior knowledge from pulsar timing,spectroscopy and distance are key to determine with certainty the most accurate re-distribution law. We propose an extension of the direct heating framework that allowsfor exploring a variety of heat redistribution effects. Future work is necessary to de-termine the relevant laws from first principles and empirically using complementaryobservations.
Key words: pulsars: individual: PSR J2215+5135 – binaries: close – convection –diffusion – stars: atmospheres
Spider pulsars are binary systems in which the primary com-ponent is a millisecond pulsar and the secondary a low-massstar, which we will call the companion in this paper. Theorbital period of the binary is typically of a few hours. Thecompanion is generally close to filling its Roche lobe andusually assumed to be tidally locked onto the neutron star.Spiders are found in two sub-species: redbacks with com-panion mass (cid:38) . M (cid:12) , and black widows with a companionmass of a few . M (cid:12) . These names were coined after twoarachnid species which share the characteristic that the lightmale companion is sometimes eaten by the heavier female. (cid:63) E-mail: [email protected];[email protected]
For their stellar counterparts, there is indeed suspicion thatthe low-mass companion is being gradually evaporated bythe intense wind of high-energy particles radiated by thepulsar (e.g. Fruchter et al. 1988). This is evidenced by thewide radio eclipses attributed to clumps of ablated mate-rial surrounding the companion far out of its Roche lobe,although it is as yet unclear whether this is sufficient to leadto the disappearance of the star (e.g. Polzin et al. 2020).This irradiation of the companion, which often exceedsthe intrinsic luminosity of the star, results in a characteris-tic day-night pattern in its light curve as it moves aroundits orbit with the pulsar. Once modelled and combined withpulsar timing and potentially with spectroscopy, optical lightcurves allow one to infer the inclination and the mass ratioof a system, and thereby the mass of the two components. Inparticular, there is evidence that the mass of spider pulsars © a r X i v : . [ a s t r o - ph . H E ] J un G. Voisin et al. could be on average larger than for other pulsars (Linares2019; Strader et al. 2019), and so they could be used to con-strain the maximum mass and equation of state of neutronstars (e.g. ¨Ozel & Freire 2016). Modelling of optical observa-tions of spider companions also provides an indirect probe ofthe pulsar wind. In particular, comparing the temperaturesof the day and night sides of the star provides an estimateof the amount of irradiating power necessary to sustain sucha difference. Modelling the interaction of the wind with thestellar material, and in particular determining what compo-nents (gamma rays, leptons, hadrons) can penetrate belowthe photosphere and produce the observable effective tem-perature difference also provides insight in the compositionof the wind (Zilles et al. 2019).Both the determinations of the orbital and wind pa-rameters are highly dependent on the modelling of the tem-perature at the surface of the companion star. A commonapproach consists in assuming a direct heating of the sur-face whereby the energy deposited by the pulsar wind isre-radiated by the companion at the exact location where itis absorbed (e.g. Breton et al. 2012). Although this approachpermits reasonable fits of some light curves (e.g. Breton et al.2013), this model is unable to account for asymmetries be-tween the leading and trailing edge of the companion (as-suming a symmetric irradiation pattern) such as seen in, forexample, the black widows PSR J2051-0827 (Stappers et al.2001) and PSR B1957+20 (Kandel & Romani 2020) or theredback PSR J2215+5135 (Romani & Sanchez 2016; Linareset al. 2018; Schroeder & Halpern 2014b).Various models have been proposed to explain asymme-tries. Two of them, the magnetic-field ducting of the ener-getic charged particles of the pulsar wind by the magneticfield of the companion (Sanchez & Romani 2017) and radia-tion from an intra-binary shock formed between the winds ofthe two components (Romani & Sanchez 2016), still assumedirect heating of the companion but change the irradiationpattern from an isotropic point source through the inter-action of the pulsar wind with the companion’s own windand/or magnetic field. A third approach consists in empiri-cally adding hot or cold spots at the surface of the star (e.g.Shahbaz et al. 2017). Thus, in all cases heat is assumed to beneither diffused nor convected within the star after energyreaches the surface.In Sect. 2 we develop a simple heat redistribution modelwithin the outer layers of the star which constitutes a nat-ural extension of direct heating models. In principle, theheat flux must be calculated from detailed stellar and atmo-spheric models of the star. However, we focus in the presentwork on demonstrating the basic properties and the inter-est of this new framework by using simple diffusion-like andconvection-like laws. In Sect. 3, we apply this simple modelto the light curve of the companion of PSR J2215+5135(Linares et al. 2018) in order to empirically determine themost probable law. We then discuss the physical interpre-tation of the results in Sect. 4. At the time of finishing thispaper, a similar model was published in (Kandel & Romani2020) which appears to be a special case of the frameworkpresented here, where no diffusion effect is considered and aparticular convection law is used. For comparison purposes,we also reproduce this model in the present work.
Currently, state-of-the-art light-curve-modelling softwaressuch as
Icarus (Breton et al. 2012) rely on the approxi-mation that the power impinging on the companion star isthermalised and re-radiated at the location on the photo-sphere where it was absorbed. This leads to the followingenergy balance, σ T = σ sb T + L w , (1)where σ sb is the Stefan-Boltzmann constant, T dh is the tem-perature of the photosphere after irradiation, T b is the basetemperature without irradiation, L w is the energy flux of thepulsar wind at the photosphere.Let us note that the base temperature T b is not nec-essarily constant over the star, but can be affected by, forexample, gravitational darkening or magnetic activity (starspots). The irradiation flux L w includes the cross-section ofthe stellar surface relative to the incoming flux. Indeed, ifthe irradiating flux is L k where k is a unit vector and thenormal to the stellar surface is given by the unit vector n then the flux crossing the surface element is L w = L k · n .The function L can take different forms depending of whatthe source of irradiation is assumed to be. It is common toassume symmetric direct heating from a point source, thatis irradiation by a wind radially expanding from the pulsar,but it has been proposed that the wind might be repro-cessed by an intra-binary shock (Romani & Sanchez 2016)or channelled by the companion’s magnetic field (Sanchez& Romani 2017) thus making L highly non-trivial in thosecases.There are several examples in the literature that showthat the direct-heating model works well when fitting someoptical light curves (e.g. van Kerkwijk et al. 2011; Bretonet al. 2013). This tells us that i) the irradiating flux is, atleast partly, deposited below the photosphere of the star asotherwise optical light curves would not be affected, and ii)that the deposition depth is probably shallow as otherwiseheat would not emerge at the entry point on the photo-sphere. These two points have recently been backed in Zilleset al. (2019) who showed that only high-energy particles( (cid:38) MeV) can deposit their energy below the photosphere,and do so at very shallow depths, typically after crossing acolumn density < / cm . In the following, we propose to supplement equation (1) byadding the possibility of energy transport within a thin shelllocated just below the photosphere. The base of the shell isassumed to be unaltered by irradiation which implies thatit is much deeper than the reach of high-energy particlesbombarding the star. The thickness of the shell must alsobe very small compared to the size of the star, and we willconsequently consider it negligible.Within this shell, we consider the stationary equationof conservation of energy, ∇ · j = (cid:219) e (2) MNRAS , 1–11 (2020) pider heat redistribution where j is the flux of energy per unit surface, and (cid:219) e is an ex-ternal source of power per unit volume which, in the presentcase, is the irradiating power of the pulsar wind. Since thisequation is linear in j , we may consider the base flux b cor-responding to the homogeneous solution, (cid:219) e = , and a par-ticular solution i corresponding to irradiation such that j = b + i . (3)The homogeneous solution b is in principle part of a generalsolution of the full set of stellar-structure equations. Makingthe simplification that the star has a spherical photosphericsurface of radius R ∗ , one has b r ( R ∗ ) = σ sb T , (4)where b r is the radial component of b and T b is the base tem-perature, that is the photospheric temperature in absence ofirradiation. This boundary condition is in fact all we needfrom the base solution for the following derivations.In order to compute the particular solution i , we use thecondition that the inner surface of the shell is unaffected byirradiation, which gives the boundary condition i r ( R i ) = , (5)where R i is the radius at the base of the shell.We now proceed to average equation (2) over the thick-ness of the shell. We start with integrating (2) over the vol-ume of an element of shell corresponding to a surface δ S atthe surface of the star between R ∗ and R i . We immediatelyobtain ∫ d V (cid:219) e = L w δ S while we can apply Gauss’ theorem tothe divergence term such that ∫ i · d S = ∫ R ∗ R i d r ∮ d C · i (cid:107) + ( i r ( R ∗ ) − i r ( R i )) δ S + (cid:13) (cid:18) ∆ RR ∗ (cid:19) , (6)where we have used the fact that R ∗ − R i = ∆ R (cid:28) R ∗ to inte-grate over a shell element of constant section δ S , d C is a lineelement of the contour of that section, and we have decom-posed the energy flux into its radial and angular components i = ( i r , i (cid:107) ) .Now, we may revert Gauss’ theorem in each 2-dimensional slice of the first term of the right-hand side ofequation (6), ∮ d C · i (cid:107) = ∫ d S r ∇ (cid:107) · i (cid:107) + (cid:13) (cid:18) ∆ RR ∗ (cid:19) , (7)where r − ∇ (cid:107) is the angular part of divergence operator inspherical coordinates, ∇ (cid:107) = θ (cid:18) ∂ sin θ∂θ u θ + ∂∂φ u φ (cid:19) , (8)where ( θ, φ ) are respectively the colatitude and longitude atthe surface of the star, and ( u θ , u φ ) are the associated unitvectors.Inserting equation (7) back into equation (6) using theboundary condition of equation (5), and differentiating withrespect to the surface elements δ S we obtain the averagedenergy conservation equation, ∇ (cid:107) · ∫ d r r i (cid:107) = − i r ( R ∗ ) + L w , (9)where, in addition, we have used the fact that ∇ (cid:107) is inde-pendent of r to take it out of the integral on the left-handside. Defining the “average” parallel energy flux as J (cid:107) = ∫ R ∗ R i d r r i (cid:107) ( r ) , (10)and introducing the irradiation temperature σ sb T = i r ( R ∗ ) ,we may rewrite equation (9) as ∇ (cid:107) · J (cid:107) = − σ T + L w . (11)The flux that escapes the star is given by j r ( R ∗ ) = b r ( R ∗ ) + i r ( R ∗ ) . Since we have made the assumption that theirradiating power is thermalised before being re-radiated,this means that the escaped flux corresponds to a black-body at temperature T ∗ such that T ∗ = T + T , (12)and that this temperature corresponds to the actual tem-perature of the plasma at the photosphere.This allows us to write our final superficial energy trans-port equation, ∇ (cid:107) · J (cid:107) = − (cid:16) σ sb (cid:16) T ∗ − T (cid:17) − L w (cid:17) . (13)One notes that if parallel energy transport can be ne-glected, that is J (cid:107) = , one naturally recovers the commondirect heating of the companion star by the pulsar wind. Inthis case, T b is directly equal to the night-side temperatureof the star. Note that here we define the night-side tempera-ture as the temperature at the point on the surface oppositeto the pulsar’s direction. This quantity is different from theeffective temperature inferred at the inferior conjunction ofthe companion, which is an average over the visible surfaceat this particular phase. We now consider that parallel energy transport follows a lawof the form J (cid:107) = − κ ∇ (cid:107) T ∗ − T ∗ f ( θ ) sin θ u φ , (14)where the first term on the right-hand side accounts fordiffusion-like effects and the second term for convection-likeeffects. The spherical coordinates are defined as for Eq. (8)with the polar axis taken to be the spin axis of the star,and the prime meridian, φ = , intersects the binary axis onthe night side of the star. The parameter κ is the diffusioncoefficient with a dimension of energy per unit temperatureper unit surface per unit time.In the convection term, we consider that the surfacetemperature T ∗ is convected by a velocity field that rotatesaround the angular-momentum axis of the star such that ifthe function f is a constant then the convecting flow is insolid rotation with a velocity field f sin θ u φ . However, thepolar convection profile f ( θ ) may be prescribed to reflecttheoretical predictions such as, for instance, equatorial jets(see, e.g., Showman & Polvani (2011) and below).Here, we have assumed that the surface temperature T ∗ is a good proxy for the transport properties of the shell. In-deed, if the energy is deposited at a shallow depth below thesurface, then parallel temperature gradients should be max-imum near the surface and so should be diffusion. Similarly,assuming a sufficiently smooth radial temperature profile inthe shell then the photospheric temperature can be chosen MNRAS , 1–11 (2020)
G. Voisin et al. as representative of convective transport. Nevertheless, thelaw of equation (14) should be considered as an effective de-scription of the physics taking place in the outer shell of thestar and not a as law derived from first principles.To go further, we assume that κ is a constant. Insertingequation (14) in the energy redistribution equation (13), weobtain κ ∇ (cid:107) T ∗ + f ( θ ) ∂ φ T ∗ = σ sb (cid:16) T ∗ − T (cid:17) − L w , (15)where ∇ (cid:107) is the angular Laplacian.Note that equation (14) is certainly not the only onepossible solution but we favour it in this article owing toits relatively mathematical simplicity while retaining someof the expected qualitative behaviour. For instance, it couldeasily be generalised to more complex convection patternsand a non-constant κ .We present in appendix A a method to solve Eq. (15).It is interesting to note that in many cases a good approx-imation can be obtained by linearising Eq. (15) around thedirect heating solution of Eq. (1) and decomposing T ∗ ontospherical harmonics. We also found that this procedure cansuccessfully be iterated in order to obtain the fully non-linearsolution to Eq. (15), thus providing a higher accuracy. Weuse the latter method in the rest of this article. The model of Eq. (15) depends on the choice of convectionprofile f ( θ ) made by the modeller based on additional theo-retical and/or empirical evidence. We have tried the follow-ing different forms, f ( θ ) = ν, (16) f ( θ ) = ν exp (cid:18) − θ w (cid:19) , (17) f ( θ ) = exp (cid:18) − θ w (cid:19) (cid:213) i = ν i H i (cid:18) θ w (cid:19) , (18) f ( θ ) = + ν if | θ | < w ; − ν otherwise . (19)In all these profiles, ν (or ν i ) is the energy flux per unittemperature transported by convection. Equation (16) cor-responds to a constant longitudinal advection, meaning thatif the properties of the superficial layer are constant acrossthe entire surface (thickness, density, thermal capacity) thena constant ν corresponds to the constant angular velocity(around the spin axis of the star) of an advection flow insolid rotation around the star. If ν > , then the flow isrotating in the same direction as the star on its orbit. Equa-tion (17) assumes that the flow is localised within a Gaussianbelt of characteristic angular width w around the equator.Equation (18) is a generalisation of Eq. (17) to an expan-sion into Hermite polynomials H n up to 3rd order. Indeed,such an expansion has been shown to be the eigen basisof the polar dependence of flow solutions to the shallow-water model developed in Showman & Polvani (2011) forsuper-rotation in atmospheres of tidally-locked exoplanets.It follows that Eq. (17) is simply Eq. (18) with ν = ν and ν i > = . Equation (19) corresponds to the particular case studied recently in Kandel & Romani (2020) if diffusion isnot included ( κ = ). In this model, a convection belt ofwidth w is rotating around the equator while matter flowswith opposite velocity at higher latitudes. All these profilesshare the property that the convection pattern is dominatedby an equatorial jet (e.g. Showman & Polvani 2011). Wenote that only Eqs. (18) and (19) include the possibility ofcounter-rotating flows. We show examples of the temperature difference with re-spect to direct heating, that is T ∗ − T dh , obtained using theabove temperature profiles of Eqs. (16)-(19) in Fig. 1. Onesees that, in every case, the changes in temperature are lo-cated near the terminator of irradiation by the pulsar, aswell as near the apex of the star in direction of the pul-sar when diffusion is enabled. This is because it is wherethe strongest temperature gradients of the direct heatingpattern are present. The additional wavy patterns that canbe distinguished are due to the limited number of sphericalharmonics used in the expansion of the solution. We havechecked that for l ≥ , these patterns entirely average outand do not bias the corresponding light curves (see nextsection).As can be seen in Fig. 1, diffusion transports energyfrom the day side to the night side symmetrically with re-spect to the binary axis (if the star is not spherical somesmall asymmetries can appear, in particular due to gravitydarkening). On the contrary, the effect of convection is asym-metric between the leading and the trailing edge of the star,and localised at particular latitudes (except for the profile ofEq. (16)). As a result, convection effectively creates hot andcold spots at the intersection of the characteristic latitude ofa stream and the terminator line. However these spots arelargely smoothed when diffusion is present. As an example, the above model was fit to the multi-colouroptical light curve ( g (cid:48) , r (cid:48) , i (cid:48) ) of the redback companion ofPSR J2215+5135 taken using the Auxiliary Port Camera(ACAM) mounted on the William Herschel Telescope in2014. This data was initially presented in Linares et al.(2018), and is also employed in Kandel & Romani (2020).This is also the same object which was used by Romani &Sanchez (2016) to demonstrate the effectiveness in invokingan intra-binary shock in order to describe the asymmetriespresent in the light curve of PSR J2215+5135.We have fitted the usual symmetric direct heatingmodel without heat redistribution, and compared it withheat redistribution models using the convection profiles ofEqs. (16)-(19) both with and without diffusion, that is with κ free or fixed to zero in Eq. (15). In each case, the param-eter space was explored using multinest (Feroz & Hobson2008; Feroz et al. 2009; Feroz et al. 2013) nested samplingalgorithm as implemented in python through pymultinest (Buchner et al. 2014). This algorithm was chosen because itallows to compare the evidence of each model, that is theprobability of the model given the data, and therefore per-form a direct comparison between them. MNRAS , 1–11 (2020) pider heat redistribution Figure 1.
Temperature maps representing direct heating corresponding to the solution of Table 1 (top-left panel) and examples oftemperature differences obtained when using the convection patterns of equations (16)-(19) without (left column) or with (right column)diffusion. The convection and diffusion parameters are chosen to serve an illustrative purpose with κ = , ν = , ν i = ν / i and w = ◦ . The blue lines show the location of the terminator line of the corresponding direct heating. Two cycles of longitudes are shownfor clarity. The point at longitude ◦ , co-latitude ◦ faces the pulsar. Apart from the heat redistribution parameters, we fitfor extinction E ( g − r ) , the amplitude of the projected radialvelocity K , distance d , base and irradiation temperatures T b and T ir , Roche-lobe filling factor f RL and system inclination i (e.g. Breton et al. 2012, 2013). We also report in Table 1some derived parameters of interest: the mass ratio q , thepulsar and companion masses M psr and M c , and the irradi-ation efficiency (cid:15) . The latter is defined as the ratio betweenthe pulsar spin-down power and the irradiating power ab-sorbed by the star (e.g. Breton et al. 2013). To derive some ofthese parameters, we made use of the orbital characteristicsobtained from pulsar timing, in particular the pulsar pro-jected semi-major axis a p sin i = . ± . lt-s andthe orbital period P = . ± . d (Abdoet al. 2013). There are three parameters for which we set informed priorswhen exploring the parameter space with
Multinest : thedistance to the source, the optical extinction in the directiontowards the source, and the radial velocity of the companionstar. In addition, the inclination had a sin ( i ) prior applied,reflecting a uniform prior on the orientation of the system.The distance prior has three components. The first isbased on the estimated space density and transverse veloc-ity of millisecond pulsars along the line of sight towardsPSR J2215+5135, with the underlying spacial density for MSPs coming from Levin et al. (2013). This component hasa Gaussian distribution in distance from the Galactic centrewith width σ = . kpc, a decaying exponential in heightabove the Galactic plane with a scale height of 0.5 kpc, anda decaying exponential in transverse velocity, with a meanvelocity of 100 km s − (Manchester et al. 2005). The sec-ond component comes from an upper limit on the system’sparallax of < . milliarcseconds at the 5 σ level, which wasobtained from the second data release of the Gaia space-craft (Gaia Collaboration et al. 2018). The third componentcomes from combining the most recent galactic electron den-sity distribution model (Yao et al. 2017) with the dispersionmeasurement value of . ± . pc cm − obtainedfrom radio timing of PSR J2215+5135. The resulting prioris shown over the relevant parameter space in the distanceplot of Fig. 4, and is not very constraining.The prior on the optical extinction, that is E ( g − r ) ,was a Gaussian centred on 0.13 and with width σ = . ,and comes from the measured value from the Bayestar19 dust maps (Green et al. 2019). The radial velocity of thecompanion star, K , had a Gaussian prior centred on 412km s − with width σ = km s − , inline with the estimatedcentre-of-mass velocity of the secondary given by Linareset al. (2018). MNRAS , 1–11 (2020)
G. Voisin et al.
It appears that models with the convection profiles of Eqs.(17)-(19) all converge to the profile of Eq. (16). Indeed, theircharacteristic width is compatible with w ∼ π . As a conse-quence, we report in detail only the results for the uniformconvection model of Eq. (16) with and without diffusion.We make an exception for the bizone convection profile of(19) without diffusion in order to compare with the recentwork of Kandel & Romani (2020). We also report for com-parison the symmetric direct heating model. The results ofthese four fits are collated in Table 1 by order of increasingevidence. The best-posterior light curves of these models arereported in Fig. 5.One sees that the most favoured model is the modelwith both uniform convection and diffusion, while uniformand bizone convection without diffusion yield quasi-identicalsolutions. According to their respective evidence, the uni-form+diffusion model is ∼ times more likely than uni-form or bizone convection alone. However, one can see thatthe ranking in terms of best χ is quite different, reflectingthe role of the priors in the results. As can be seen in Fig. 4,the uniform convection+diffusion model is the only one thatfits the best within the distance and extinction priors, com-parably to the direct heating model, while the two purelyconvective models stand at the edge of the extinction priorand require a significantly larger distance.This correlates with the fact that these models requireboth very high base temperature T b (cid:39) K and irradiationtemperature T ir (cid:39) K implying a maximum day-side tem-perature over T ( max ) D ∼ ( T + T ) / (cid:39) K. Spectroscopicobservations reported by Linares et al. (2018) provide av-erage night and day-side temperatures of T ( spec ) N = + − and T ( spec ) D = + − K respectively. These temperatures arederived from spectra taken at inferior and superior conjunc-tion of the companion respectively. These spectra result fromthe superposition of light originating from within the visiblesurface of the star which is not at a uniform temperatureand therefore should be seen as average values. In particu-lar, they are not equal to the minimum night-side tempera-ture ( ∼ T b ) and maximum day-side temperature T ( max ) D . Wehave estimated T ( spec ) N and T ( spec ) D for our models by com-puting the position of the peak of the spectrum resultingfrom the sum of the local black-body spectra of each visiblesurface elements at inferior and superior conjunction respec-tively. The results, reported in Table 1 show that only theuniform convection+diffusion model, and with slightly moretension the symmetric direct heating model, are compatiblewith the spectroscopic observations of Linares et al. (2018),while the convection-only models require much larger tem-peratures for both sides of the star. We note the very impor-tant role of diffusion here. Indeed, diffusion simultaneouslydecreases the day-side temperature and increases the night-side temperature by ∼ K compared to direct heatingwith the same parameters, as is shown on Fig. 2. It entailsthe much smoother temperature map of Fig. 3 compared to,for instance, the direct heating model (top left panel of Fig.5) which allows a moderate day-night temperature differencedespite the significantly larger irradiation temperature andcooler base temperature.
Figure 2.
Maps of temperature difference with respect to directheating for the best posterior parameters of the bizone convec-tion model (top), the uniform convection model without diffusion(middle) and with diffusion (bottom) from the results of the light-curve fits of PSR J2215+5135’s companion presented in Table 1.The blue lines show the location of the terminator line of the cor-responding direct heating pattern. Two cycles of longitudes areshown for clarity. The point at longitude ◦ , colatitude ◦ facesthe pulsar. Figure 3.
Temperature map of the best-posterior solution usingwith the uniform convection+diffusion model from the results ofthe light-curve fits of PSR J2215+5135 companion presented inTable 1. The blue line shows the location of the terminator line ofthe corresponding direct heating pattern. Two cycles of longitudesare shown for clarity. The point at longitude ◦ , colatitude ◦ faces the pulsar. Interestingly, the base temperature of the uniform convec-tion+diffusion model is significantly lower than any of theother models with T b = + − K in comparison to at least
K for the symmetric direct heating model and morethan
K for the convection-only models. Let us remem-ber that the base temperature is the effective temperaturethat the star would have in absence of irradiation. Due to thelack of knowledge of the stellar structure of redback compan-ions, it is difficult to know what this temperature should bein theory. The net result is that, especially for cases wherethe heat is substantially redistributed to the back of thestar, the observed average night side temperature may de-part significantly from the true temperature that it wouldhave without irradiation.
MNRAS , 1–11 (2020) pider heat redistribution Direct Bizone Uniform Uniform + Diffusion log Z − − . − . − . N dof
229 227 228 227 χ . . . χ . . . Fitted parameters E ( g − r ) . + . − . . + . − . . + . − . . + . − . K (km/s) + − + − + − + − d (kpc) . + . − . . + . − . . + . − . . + . − . T b (K) + − + − + − + − T ir (K) + − + − + − + − f RL . + . − . . + . − . . + . − . . + . − . i ( ◦ ) . + . − . . + . − . . + . − . . . + . − . κ (W/K/m ) - - - + − ν (W/K/m ) - + − + − + − w (rad) - . + . − . - -Derived parameters q . + . − . . + . − . . + . − . . + . − . M psr ( M (cid:12) ) . + . − . . + . − . . + . − . . + . − . M c ( M (cid:12) ) . + . − . . + . − . . + . − . . + . − . (cid:15) . + . − . . + . − . . + . − . . + . − . T ( spec ) N (K) + − + − + − + − T ( spec ) D (K) + − + − + − + − Table 1.
Evidence sampling results for the three main models applied to J2215+5135: Direct heating, uniform convection withoutdiffusion, and uniform convection with diffusion. log Z is the natural logarithm of the model evidence, and models are ranked by increasingevidence. N dof is the number of degrees of freedom of each model, and we give the χ of the solution with the best likelihood (but notnecessarily the best posterior probability) and of the median solution. Model parameters are reported for the median solution with the95% confidence interval ( ± . ). It is however interesting to compare redback compan-ions to companions of cataclysmic variables. Indeed, thesestars have similar masses to redback companions, and sim-ilarly underwent Roche-lobe filling and mass transfer tothe benefit of their primary (a white dwarf in this case).Nonetheless, these stars are not irradiated, letting their basetemperature being seen, and have been largely studied bothobservationally and theoretically. One may therefore spec-ulate that their effective temperature is similar to the basetemperature of redback companions, in which case it wouldappear to be in the range − K depending on themass of the star (Knigge et al. 2011). Due to its large un-certainty on the lower bound of the base temperature, theuniform convection+diffusion model is the only model com-patible with this range.
The framework proposed in section 2 allows to redistributeenergy at the surface of the star assuming a given trans-port law. However, determining such a law from first princi-ple requires to determine not only the relevant microphysicsbut the hydrodynamical properties of stellar matter as well.This is out of the scope of the present paper and shouldbe addressed in future work. Here we have focused on thestudy of the simplest possible convection-like and diffusion-like redistribution laws, by assuming only a latitudinal de-pendence for the convection profile (see Sect. 3.1) and a constant diffusion coefficient. In the following we derive or-ders of magnitude to show that the values we obtain forthe convection and diffusion parameters of the uniform con-vection+diffusion model applied to PSR J2215+5135’s com-panion, ν and κ respectively, can be compatible with somesimple physical processes. The fact that all convection profiles converged to a solutionsimilar to the uniform rotation profile may simply mean thatthe finer details cannot be resolved with the available data.In particular, any latitudinal structure is bound to be largelyif not completely averaged out since photometric informa-tion only provides the total flux contribution as a functionof rotational phase with very little handle on the other axis.One may tentatively interpret the value of the convectionparameter in terms of a wind velocity similar to the one-dimensional model of Cowan & Agol (2011). Thus, assuminguniform rotation of a shell of uniform column density Σ atangular velocity ω one can write ν = Σ c p ω, (20)where c p = k b / µ (cid:39) − kg − is the specific heatcapacity of a perfect gas at constant pressure, k b is Boltz-mann’s constant, and we have approximated the meanmolecular mass µ to the mass of a proton. We may adoptthe median value of the depth of maximum heat deposition, / cm (see Zilles et al. (2019) and section 2), as a fiducial MNRAS , 1–11 (2020)
G. Voisin et al. . . . . . d ( k p c ) T b ( K ) T i r ( K ) . . . . f R L i () ( W / K / m ) ( W / K / m ) K ( m s ) . . . . . M p s r ( M ) . . . . M c ( M ) . . . . q .
04 0 .
08 0 .
12 0 .
16 0 . E ( g r ) . . . . .
50 2 . . . . . d (kpc) T b (K) T ir (K) .
75 0 .
78 0 .
81 0 . f RL
64 72 80 88 i ( ) (W/K/m ) (W/K/m )
400 408 416 424 K (m s ) . . . . . M psr ( M ) .
24 0 .
28 0 .
32 0 . M c ( M ) .
75 6 .
90 7 .
05 7 . q .
50 0 .
75 1 .
00 1 .
25 1 . Figure 4.
Results of the multinest fit of the light curve of PSR J2215+5135 using the uniform convection+diffusion model. The solidlines in the plots along the diagonal show the prior functions used, if a prior was specified. value for Σ . This assumes that deeper layers are not affectedby latitudinal convection. The characteristic hydrodynam-ical velocity is the speed of sound c s ∼ (cid:112) k b T ∗ / µ giving afiducial ω = c s / R ∗ . It follows that ν = − m − (cid:18) Σ
500 g / cm (cid:19) (cid:18) T ∗ (cid:19) / (cid:18) . R (cid:12) R ∗ (cid:19) , (21) where we have derived the stellar radius R ∗ from the resultof the fit. Remarkably, the above fiducial value for ν agreesin order of magnitude with the results of Table 1.Interestingly, the bizone model does not reproduce thefit reported recently in Kandel & Romani (2020) who usesthe same model and the same photometric data. In partic-ular, in Kandel & Romani (2020) it is found that w (cid:39) ◦ ( θ c in their notations) while we find w > ◦ , renderingthe model virtually equivalent to uniform convection. It isunclear why this happens, however we note that the othermajor difference in the results of Kandel & Romani (2020) MNRAS , 1–11 (2020) pider heat redistribution M a g n i t u d e Direct Heating R e s i d ( m a g ) Convection M a g n i t u d e Convection + Diffusion R e s i d ( m a g ) Figure 5.
The light curve of PSR J2215+5135 in SDSS g (cid:48) (green), r (cid:48) (red), and i (cid:48) (orange) bands (top panel) and residuals aftersubtraction of the model light curve (bottom panel). In each case, the best posterior models is shown. is the relatively mild base temperature T b ( T N in their no-tations) as well as irradiation temperature. These quantitiescorrelate with the level of extinction, and indeed one cansee that the fitted values of E ( g − r ) for our bizone and uni-form models (Table 1) are substantially into the tail of ourprior on this parameter (see Sect. 3.3.1). On the other handthe fit of Kandel & Romani (2020) assumes a fixed value ofextinction (corresponding to the centre of our prior). Otherreasons for the discrepancy might include their addition of aveiling flux (which we cannot assess with only photometricdata) although they report that the inclusion of this extracomponent improves the fit without affecting the fitted pa-rameters. In the outer stellar envelope the main diffusion mechanism isradiative diffusion whose flux is j rad = −( / ) σ sb T l p ∇ T (e.g.Kippenhahn et al. 2012), where l p is the photon mean freepath in the material and σ sb is Stefan-Boltzmann’s constant.Taking the average defined in equation (10) of the radiativediffusion flux over a slab of stellar matter of height H (cid:28) R ∗ ,we can estimate J (cid:107) rad ∼ − κ rad ∇ (cid:107) T ∗ , (22)where κ rad ∼ σ sb T ∗ l p HR ∗ . (23)The photon mean free path depends on the complex inter-play of density, temperature and molecular composition (e.g.Kippenhahn et al. 2012). Consequently, it is not possible tohave a precise estimate of l p without a full modelling of atleast the outer layer of the star.We note that l p = ( k ρ − ) where ρ is the local densityand k the opacity of the material. In addition, H ∼ Σ / ρ where, as before, Σ is the corresponding column density ofthe slab. Inserting a typical value for the opacity in Eq. (22), that is k ∼ / g , we can estimate the density nearthe photosphere necessary to obtain a given value of thediffusion coefficient κ rad , ρ ∼ × − g / cm (cid:18) κ rad × W K − m − (cid:19) − / (24) (cid:18) k / g (cid:19) / (cid:18) Σ
500 g / cm (cid:19) / (cid:18) R ∗ . R (cid:12) (cid:19) − (cid:18) T ∗ (cid:19) / . This value is an order of magnitude smaller than the pho-tospheric density of Solar-type stars (e.g. VandenBerg et al.2008) which have a similar surface gravity log g (cid:39) . andtemperature. However, this estimate of surface gravity doesnot take into account the effect of the star being close to fill-ing its Roche lobe, which is bound to diminish the effectivegravity near the surface of the star and in the atmospherecompared to the isolated case, thus diminishing the pres-sure and the density at the photosphere. We also note thehigher temperature (compared to the Sun) on the day sideof the companion, which would also tend to decrease thedensity at equal pressure. Thus, it seems possible that thevalue of the diffusion coefficient resulting from our fit canbe explained with radiative diffusion in the outer layers ofthe star, although a complete modelling of its atmosphericand sub-photospheric structure is necessary to answer thisquestion with certainty. Pulsar timing of PSR J2215+5135 measures the projectedsemi-major axis of the pulsar as well as its orbital period.Using Kepler’s third law, one combines these two parame-ters to compute the value of the so-called mass function (e.g.Lyne & Graham-Smith 2012) which relates the two masses ofthe system to the orbital inclination. In order to lift the de-generacy between masses and inclination, one needs two ad-ditional measurements. In the case of PSR J2215+5135, themass ratio can be inferred from spectroscopic measurements
MNRAS , 1–11 (2020) G. Voisin et al. of the companion’s projected radial velocity amplitude K ,though with extra complications due to the irradiation ef-fects (Linares et al. 2018). On the other hand, fitting theoptical light curve allows one to measure the inclination.Our results in Table 1 show that this quantity is highlymodel-dependent, ranging from . + . − . ◦ for the directheating model to . + . − . ◦ for the diffusion+uniform convec-tion model through somewhat intermediate values for thetwo convection-only models. Accordingly, the pulsar massranges from . + . − . to . + . − . for the direct heating anddiffusion+uniform convection respectively. The direct heat-ing values confirm those found in Linares et al. (2018) usingthe same dataset, while being hardly compatible at the 95%level with the diffusion+uniform convection values. Interest-ingly, the latter gives a similarly high inclination to what wasfound in Romani et al. (2015) (see Linares et al. (2018) fora review of previous measurements). In Kandel & Romani(2020), it is however argued that this previous result mighthave been biased by an extra blue veiling flux at the epochof the observations of Schroeder & Halpern (2014a) (whoseoptical light curve they use), as suggested by a correspond-ing excessive night-side temperature of that fit compared tothe spectroscopic constraints of Linares et al. (2018).In the present work, we see that inclination is substan-tially changing from one model to another, everything elsebeing equal. Although we cannot here assess with certaintythe presence of a veiling flux for lack of spectroscopic obser-vations, the night-side temperature of the diffusion+uniformconvection fit is not excessively large as discussed in Sect. 3.Another possible caveat is the lack of a significant portion ofthe orbital light curve (see Fig. 5), which might bias the fitespecially considering the asymmetry of the light curve. Wetherefore conclude that a thorough investigation involvingsimultaneous spectroscopy and photometry across an entireorbit is desirable in order to be able to reduce the risk ofbias. In this paper we have considered the effects of heat redistri-bution at the surface of companion stars of spider pulsars. Ineffect, we have supplemented the usual direct heating modelof irradiation, Eq. (1), with a single extra term accountingfor the divergence of the heat flux within the stellar surface,Eq. (13). This may be seen as the simplest addition possi-ble to direct heating models. On the other hand, the heatflux itself requires a complex modelling of the outer layerand atmosphere of the star combining microphysics, hydro-dynamics and thermodynamics which is outside of the scopeof the present work.In the spirit of studying the simplest possible exten-sions to direct heating models we have evaluated the effect ofsimple convection-like and diffusion-like laws, Eqs. (14) and(15). The solution of the redistribution equation can be rep-resented under the form of heat redistribution temperaturemaps, Fig. 1, which show that both convection and diffu-sion effects are most intense near the irradiation terminatoror near the apex of the star towards the pulsar. Interest-ingly, convection is naturally able to produce patterns akinto hot or cold spots at the terminator. We also note thatheat redistribution models are compatible with other mod- els which modify the irradiation pattern such as intra-binaryshock models (Romani & Sanchez 2016) or magnetic-fieldducting models (Sanchez & Romani 2017), and with modelsthat modify the base temperature such as hot and cold spots(e.g. Shahbaz et al. 2017).We have applied our models to the light curve of the al-ready well-studied companion of the redback pulsar PSRJ2215+5135 in order to determine empirically the mostlikely form of the heat flux. Various convective flows withand without diffusion were tried. We found that, althoughevery redistribution model provides a substantially better fitthan the symmetric direct heating model, the model asso-ciating diffusion to convective flows in uniform rotation ismost likely (see Table 1).However, since with every model substantial fit residu-als remain these results should be taken with caution and weconsider that the main value of the different fits lies in thecomparison with each other. Indeed, as it appears in Table 1,the various models lead to sometimes very discrepant fittedparameters, in particular concerning the base and irradia-tion temperatures, the inclination, the filling factor or theirradiation efficiency. This suggests that, on top of detailedmodelling, the determination of the “true” model of heatredistribution will certainly require complementary obser-vations such as spectroscopic measurements of the effectivetemperature, or accurate and independent distance measure-ments.
ACKNOWLEDGEMENTS
The authors acknowledge support of the European ResearchCouncil, under the European Union’s Horizon 2020 researchand innovation program (grant agreement No. 715051; Spi-ders).
REFERENCES
Abdo A. A., et al., 2013, ApJS, 208, 17Breton R. P., Rappaport S. A., van Kerkwijk M. H., Carter J. A.,2012, The Astrophysical Journal, 748, 115Breton R. P., et al., 2013, The Astrophysical Journal, 769, 108Buchner J., et al., 2014, A&A, 564, A125Cowan N. B., Agol E., 2011, The Astrophysical Journal, 726, 82Feroz F., Hobson M. P., 2008, MNRAS, 384, 449Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2013, arXive-prints, p. arXiv:1306.2144Fruchter A. S., Stinebring D. R., Taylor J. H., 1988, Nature, 333,237Gaia Collaboration et al., 2018, A&A, 616, A1Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D.,2019, ApJ, 887, 93Kandel D., Romani R. W., 2020, arXiv:2002.12483 [astro-ph]Kippenhahn R., Weigert A., Weiss A., 2012, Stellar Struc-ture and Evolution. Astronomy and Astrophysics Library,Springer Berlin Heidelberg, Berlin, Heidelberg, http://link.springer.com/10.1007/978-3-642-30304-3
Knigge C., Baraffe I., Patterson J., 2011, The Astrophysical Jour-nal Supplement Series, 194, 28Levin L., et al., 2013, MNRAS, 434, 1387Linares M., 2019, arXiv:1910.09572 [astro-ph]Linares M., Shahbaz T., Casares J., 2018, The Astrophysical Jour-nal, 859, 54 MNRAS , 1–11 (2020) pider heat redistribution Lyne A., Graham-Smith F., 2012, Pulsar Astronomy, 4 editionedn. Cambridge University Press, Cambridge ; New YorkManchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ,129, 1993Olver F. W. J., National Institute of Standards and Technology(U.S.) eds, 2010, NIST handbook of mathematical functions.Cambridge University Press : NIST, Cambridge ; New York¨Ozel F., Freire P., 2016, Annual Review of Astronomy and Astro-physics, 54, 401Polzin E. J., Breton R. P., Bhattacharyya B., Scholte D., SobeyC., Stappers B. W., 2020, Monthly Notices of the Royal As-tronomical SocietyRomani R. W., Sanchez N., 2016, ApJ, 828, 7Romani R. W., Graham M. L., Filippenko A. V., Kerr M., 2015,The Astrophysical Journal, 809, L10Sanchez N., Romani R. W., 2017, ApJ, 845, 42Schroeder J., Halpern J., 2014a, The Astrophysical Journal, 793,78Schroeder J., Halpern J., 2014b, ApJ, 793, 78Shahbaz T., Linares M., Breton R. P., 2017, Monthly Notices ofthe Royal Astronomical Society, 472, 4287Showman A. P., Polvani L. M., 2011, The Astrophysical Journal,738, 71Stappers B. W., van Kerkwijk M. H., Bell J. F., Kulkarni S. R.,2001, The Astrophysical Journal, 548, L183Strader J., et al., 2019, The Astrophysical Journal, 872, 42VandenBerg D. A., Edvardsson B., Eriksson K., Gustafsson B.,2008, The Astrophysical Journal, 675, 746Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29Zilles A., Kotera K., Rohrmann R., Althaus L., 2019, MonthlyNotices of the Royal Astronomical Societyvan Kerkwijk M. H., Breton R. P., Kulkarni S. R., 2011, TheAstrophysical Journal, 728, 95
APPENDIX A: SOLUTION OF THE HEATREDISTRIBUTION EQUATION
In this section we use the lighter notation T ≡ T ∗ . A1 Solution of the linearised transport equation
Equation (13) is strongly non-linear in T due to the T term.However, given the relative success of direct heating modelswe may assume that energy redistribution is only a pertur-bation of the temperature distribution at the surface of thestar. Thus, we may write T = T dh + t and, assuming t (cid:28) T dh ,expand equation (13) to first order in t , at − κ ∇ (cid:107) t − f ( θ ) ∂ φ t = s , (A1)where a = σ sb T , (A2) s = κ ∇ (cid:107) T dh + f ( θ ) ∂ φ T dh . (A3)The linearised equation (A1) can be solved algebraicallyafter decomposing the function onto the orthogonal basis ofspherical harmonics { Y lm } l ≥ − l ≤ m ≤ l (e.g. Olver & NationalInstitute of Standards and Technology (U.S.) 2010). In thisbasis, the functions t , a and s are represented by the vectors t , a and s respectively, such that t = (cid:213) t lm Y lm , (A4) where t lm are the coefficients of t , and similarly for s , a and f . It follows that equation (A1) can be expanded into (cid:213) lml (cid:48) m (cid:48) t lm a l (cid:48) m (cid:48) Y lm Y l (cid:48) m (cid:48) + κ (cid:213) lm l ( l + ) t lm Y lm − (cid:213) lml (cid:48) m (cid:48) imt tm f l (cid:48) m (cid:48) Y lm Y l (cid:48) m (cid:48) = κ (cid:213) lm s lm Y lm . (A5)By projecting equation (A5) onto each spherical harmonicwe obtain a set of linear algebraic equations the solution ofwhich is formally given by t = M − s , (A6)where we have introduced the matrix M = { M ij } . Its coeffi-cients are defined by M αβ = (cid:213) γ µ αβγ ( a γ − im β f γ ) + κ l α ( l α + ) δ αβ , (A7)where δ αβ = if α = β and otherwise, and eachindex α, β, γ maps onto a different pair of spherical-harmonic indices ( l , m ) (for example α = { , , , ... } →( l α , m α ) = {( , ) , ( , − ) , ( , ) , ( , ) ... } ). We have introducedthe spherical-harmonic multiplication coefficients { µ αβγ } such that Y β Y γ = (cid:213) α Y α µ α,β,γ . (A8)These coefficients can, for example, be obtained from theClebsch-Gordan coefficients (e.g. Olver & National Instituteof Standards and Technology (U.S.) 2010). A2 Solution of the full non-linear transportequation
Some irradiated stars show very large temperature differ-ences between their day and night sides, to the point thatthe temperature difference might exceed the temperature ofthe night side. In this case, the assumption that heat redis-tribution is only a perturbation of direct heating may fail.Here, we propose a fixed-point scheme to solve the full non-linear equation (15) by iterating the linearised solution ofsection A1.At each iteration, the temperature distribution T n + iscalculated according to T n + = T n + t n + , (A9)where t n + is the solution of equation (15) linearised withrespect to T n , A n t n + − κ ∇ (cid:107) t n + − f ( θ ) ∂ φ t n + = S n , (A10)where A n = σ T n , (A11) S n = κ ∇ (cid:107) T n + f ( θ ) ∂ φ T n − σ ( T n − T ) + L w . (A12)Equations (A9) and (A10) form a sequence that can beinitialised with T = T dh , t = such that t is equal to t ofthe previous section. The solution of equation (A10) is givenby equations (A6) and (A7) only replacing the vectors a , s , t by the corresponding a n , s n , t n + .In practice, this scheme converges after a few iterationswith the stopping criterion (cid:107) t n + (cid:107) < K. This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000