Mon.Not.R.Astron.Soc.000,??–??(2011) Printed1February2012 (MNLATEXstylefilev2.2) Phase dependent view of Cyclotron lines from model accretion mounds on Neutron Stars Dipanjan Mukherjee 1, Dipankar Bhattacharya 2 2 1 1 IUCAA, Post Bag 4, Pune 411007, India. Email: [email protected] 0 2 IUCAA, Post Bag 4, Pune 411007, India. Email: [email protected] 2 n a 1February2012 J 1 3 ABSTRACT In this paper we make a phase dependent study of the effect of the distortion of ] E local magnetic field due to confinement of accreted matter in X-ray pulsars on the H cyclotron spectra emitted from the hotspot . We have numerically solved the Grad- Shafranovequation for axisymmetric static MHD equilibria of matter confined at the . h polar cap of neutron stars. From our solution we model the cyclotron spectra that p will be emitted from the region, using a simple prescription and integrating over the - entire mound. Radiative transfer through the accretion column overlying the mound o may significantly modify the spectra in comparison to those presented here. However r t we ignorethis in the present paper in order to expose the effects directly attributable s a to the mound itself. We perform a spin phase dependent analysis of the spectra to [ study the effect of the viewing geometry. 2 Keywords: accretion—line:formation—magneticfields—radiationmechanisms: v non-thermal — (stars:) binaries: general — X-rays : binaries 0 5 8 2 . 1 INTRODUCTION pressure from the confined plasma will be reflected in the 0 spectra emitted from theboundary of thisregion. 1 Neutron stars in high mass X-ray binaries have high mag- 1 netic fields ( 1012G) and accrete matter from their com- The nature and variation of the cyclotron spectra can 1 panion stars∼either via stellar winds or by disc accre- giveimportantcluesregardingthepropertiesoftheemission : v tion. Magnetospheric interaction with the accretion flow region.Manysystemsshowvariationsoflineenergiesofthe i causes the matter to be channelled to the magnetic CRSFwiththephaseofrotatione.g.VelaX-1,HerX-1,4U X poles, forming accretion columns (see e.g. Ghosh et al. 0115+63,GX 301–2 etc. This can bedueto thevariation of r the local magnetic field structure at one or both poles as a (1977), Ghosh & Lamb (1978), Koldoba et al. (2002) and a line of sight moves across the neutron star. Apart from the Romanovaet al. (2003)). The infalling plasma, with ini- spinphasedependence,thecyclotronspectraarealsoseento tial relativistic infall velocities, passes through an accretion dependontheluminositystateofthesystem.Somesystems shock at a height of a few kilometres from the neutron star likeV0332+53(Tsygankov et al.2010)showanegativecor- surface and then settles down to a gradually slowing sub- relation betweenluminosityandcyclotronlineenergywhile sonic flow (Brown & Bildsten 1998; Cumming et al. 2001). somelikeHerX-1(Staubert et al.2007)showapositivecor- Such X-ray binary systems show characteristic cy- relation. Suchdependenceofthelineenergywith changein clotron resonance scattering features (CRSF) in their spec- accretionratesuggestschangeoflocalgeometryormagnetic tra resulting from resonant scattering of radiation by elec- field structure. Some sources (e.g. 4U 1538–52, A 0535+26, trons in the presence of strong magnetic field (for discus- V 0332+53 etc) show multipleabsorption features with an- sion on theory and observation of cyclotron scattering fea- harmonicseparation whichcanbeduetodistortion oflocal tures see e.g. Harding & Preece (1987), Araya& Harding field from dipolar magnetic field (Nishimura 2005, 2011). (1999), Araya-G´ochez& Harding (2000), Sch¨onherret al. (2007) and Mihara et al. (2007)). In the immediate post In this paper we examine the effect on the cyclotron shockregion theflowvelocities arestill relativistic ( 0.16c spectra arising from the distortion of local magnetic field ∼ forγ =5/3gas)andtheplasmaisopticallythintocyclotron caused by the confined plasma. We consider an accreted scattering. As the accreted plasma descends and cools, it mound in static equilibrium confined by the magnetic field forms at the base a static mound confined by the magnetic at the magnetic pole of a neutron star. We construct the field, and becomes optically thick to cyclotron scattering. equilibrium solution by solving the Grad-Shafranov equa- Any distortion of the magnetic field in the mound due to tion. We do not consider the effects of continued accretion (cid:13)c 2011RAS 2 Mukherjee and Bhattacharya in this paper. We model the X-ray emitting hotspot as a 2 STATIC MHD EQUILIBRIA OF ACCRETED mound of accreted matter with finite height and no atmo- MATTER ON NEUTRON STAR POLES sphere.TheGrad-Shafranov equationfor theaccreted mat- Inthisworkweconsideraneutronstarbinarysystemwhere ter on theneutron star poles has been previously solved by the magnetosphere cuts off the accretion disc at Alfv´en otherauthorse.g.Hameury et al.(1983),Brown & Bildsten radius (Ghosh et al. 1977; Ghosh & Lamb 1978) and mat- (1998), Litwin et al. (2001), Melatos & Phinney (2001), ter is accreted on to the polar cap. We will consider a Payne& Melatos (2004), Payne& Melatos (2007) and Vigelius & Melatos(2008)whosemainaimwastostudythe typical slowly spinning neutron star of mass 1.4M⊙, ra- dius R = 10 km and magnetic field B = 1012 G (e.g. extent of deformation and stability of the confined mound Bhattacharya & van den Heuvel(1991)).Apolar capofra- and also to deduce the effects of magnetic screening on the dius R = 1 km will be considered corresponding to the dipole moment of a neutron star. In this paper we extend p footprint of dipole field lines that extend beyond a typical thisbodyofworktopredictthecyclotronspectraemanating Alfv´enradiusof 1000 km.Accretedmatterisassumed to from such mounds. ∼ be confined within the polar cap region. We will consider the accreted matter to form a degenerate mound of finite We adopt a geometry similar to that used by height( 100morless)withapolytropicequationofstate. Hameury et al. (1983), Brown & Bildsten (1998) and ≃ Weassume that themound is in a steady state equilibrium Litwin et al.(2001)andanumericalalgorithmsimilartothe supportedbythemagneticfield.Weworkinacylindricalge- onedevelopedbyMouschovias(1974)andPayne& Melatos ometry(r,θ,z)with origin at thebase ofthepolar cap (see (2004) (PM04) for solving the Grad-Shafranov equation. Appendix B)andconsiderNewtoniangravitywithconstant HoweverourtreatmentdiffersfromPM04inseveralaspects. acceleration (Hameury et al. 1983; Litwin et al. 2001). We work in an axisymmetric cylindrical coordinate system insteadofthesphericalcoordinatesystemofPM04.Weuse g= 1.86 1014 M∗ Rs −2 cm s−2 zˆ (1) a polytropic equation of state for the accreted gas instead − × 1.4M⊙ 10km of the isothermal equation of state of PM04. Finally, we (cid:18) (cid:19)(cid:18) (cid:19) consider the mound to bestrictly confined to the polar cap The initial magnetic field (when no accreted matter is region, while PM04 allowed a significant amount of mass present) is dipolar. We approximatethe dipolar field in the loading outside thepolar cap. region by an uniform field along z : Bd = B0ˆz. We as- sumeaxisymmetryofthepolarcapmoundandusetheideal MHDequations,whichmaybecastintheformoftheGrad- Wesimulatethecyclotronspectraemittedfromtheac- Shafranov equation. We solve this numerically to find the creted mound and perform a phase resolved analysis of the fieldandmatterdensityconfiguration forthestaticequilib- emission. Our main objective in this paper is to perform a rium solution of the system. phase dependent study of the effects of accretion induced distortion of the local magnetic field on the emergent spec- tra. In this work we do not perform a detailed radiative 2.1 Formulation of Grad-Shafranov equation transfer calculation of CRSF. Instead we use a Gaussian profilefor thecyclotron feature originating from each point For an axisymmetric system, one may decompose the mag- oftheemission region,withthecentrallineenergygivenby neticfield into a poloidal and a toroidal part : the magnetic field strength at that point, according to the ∇ψ θˆ well known relativistic formula given by Sokolov & Ternov B=Bp+Bθ = r× (2) (1968).Wealsoincorporatetheeffectsofgravitationalbend- ing of light and finite energy resolution of detectors. We For our work we will assume B = 0. The function ψ(r,z) θ generate the resultant spectra by integrating over the en- is the flux function which at a fixed r and z is proportional tire mound, taking into account the variation of the field tothepoloidal fluxpassing through a circle of radius r(see strength over theemitting region. Kulsrud(2005);Biskamp(1993)formorediscussiononthis). The poloidal components of themagnetic field are We structure the paper as follows. In Sec. (2) we first 1∂ψ 1∂ψ review the formulation of the Grad-Shafranov equation for Br =−r ∂z ; Bz = r ∂r (3) thestaticMHDequilibria.Wethenoutlinethenumericalal- Using eq. (2) we can write thestatic Euler equation as gorithm adopted to solve theGrad-Shafranov equation and the test cases for verifying the code. In Sec. 3 we discuss ∆2ψ the nature of the solutions obtained by solving the Grad ∇p ρg+ ∇ψ=0 (4) − 4πr2 Shafranovequationanddiscusstherangeofparameterspace withinwhichthevalidsolutioncanbeobtained.Ourresults where∆2istheGrad-Shafranovoperator:∆2 =r∂∂r(r1∂∂r)+ areindicativeof theonset ofMHD instabilities beyondthis ∂∂z22. Assuming an adiabatic gas p = kadργ we sep- boundary.InSec. 4wedescribethealgorithmusedtosimu- arate the r and z components (Hameuryet al. 1983; latethespectrafromthemoundanddiscusstheresultsfrom Litwin et al. 2001) by the method of characteristics (sim- our simulation of the cyclotron absorption features. In Sec. ilar to Payne & Melatos 2004) 5 we summarise the results obtained from the simulations ∆2ψ ∆2ψ of the spectra and discuss the implications on observations p ρg+ ψ =0 ; p + ψ =0 (5) z− 4πr2 z r 4πr2 r of actual sources. The technical details of the geometrical constructusedtocomputethespectraarepresentedin Ap- where the subscripts indicate partial derivatives. Eliminat- pendix B. ing ∆2ψ from (5) we get the equation of the integral curve 4πr2 (cid:13)c 2011RAS,MNRAS000,??–?? Phase dependent view of Cyclotron lines 3 as : γ vs Density (in g cm−3) 1.7 dr dρ dz= = (6) −ψz/ψr −ρg/c2s ρ 1.6 n awbhoevreetcw2so=eγqup/aρtioisntsh,ewaedgieatbψati=cscpoenesdtaonfts(owunhdic.hSmolevainngstthhee np/dl 1.5 dl solutions are on constant ψ surfaces) and gz + γp = = (γ−1)ρ γ 1.4 f(ψ). Here f(ψ) is a ψ dependent constant of integra- tion. Rearranging the terms we can write the density as 1.3 (Hameury et al. 1983) 100 102 104 106 108 Density in g cm−3 1 ρ=A[Z0(ψ) z]γ−1 (7) − Figure 1. Plotofadiabatic index(γ) vsdensity (gcm−3).For 1 whereA=[g(γ 1)/(γkad)]γ−1 isaconstant.Thefunction lower densities γ asymptotically converges to 1.6667, which is − Z0(ψ)isthemoundheightprofilewhichdefinesthevertical thevalueobtainedinthenon-relativisticapproximation.Forρ≥ height of the mound for an adiabatic gas expressed in flux 108gcm−3 γ converges to 1.333 which is the value obtained in coordinate (ψ) instead of r. The values of ρ,p and their theultra-relativisticapproximation. derivativesgo to zero smoothly at z=Z (ψ). 0 Puttingeq. (7)ineq. (4)weobtaintheGrad-Shafranov AsshowninFig. 1,atρ 107g cm−3 wehaveγ 1.4.For equation (hereafter G-S) for an adiabatic gas. ≃ ≥ densities lower than this γ rises and reaches 1.667 asymp- ∆2ψ dZ totically.Since,forρ 107 g cm−3,γ issignificantlyhigher = ρg 0 (8) ≤ 4πr2 − dψ than 1.33 we model the accreted matter in the mound as a degeneratenon-relativistic electron gas with athermal pro- The Grad-Shafranov equation is a coupled non-linear ellip- ton background whose pressure is negligible compared to tic partial differential equation. We have solved the Grad- electrondegeneracypressure.Theequationofstateforsuch Shafranovequationnumericallyfollowingthealgorithmout- a system is (valuesquoted are in cgs) lined in Appendix A. ¯h2 ρ 5/3 p=[(3π2)2/3 ] =3.122 1012ρ5/3 (10) 5me (cid:18)µemp(cid:19) × 3 RESULTS AND ANALYSIS OF SOLUTIONS The plasma is dominated by degeneracy pressure if OF GRAD-SHAFRANOV EQUATION T <1, where T is theFermi temperature : TF F Wehavemadeseveralrunswithdifferentmoundheightpro- files.Thesolutionsshowexpectedbehaviourofmatterpush- T = mec2[ X2 +1 1] ing field lines outwards until the tension in the field lines F KB F − q support thegas pressure. The deformation of thefield lines X = pF = 1 ( 3h3 )1/3ρ1/3, p being the Fermi mo- increasewithlargerbasepressureanddensityofthemound F mec mec 8πµemp F mentum. Observed X-ray continuum from high mass X-ray till the solution breaks down after a threshold density (see binaries(Coburn et al.2002;Becker & Wolff2007)indicate Sec. 3.3).Inthissectionwediscussthesolutionsfromsome that photospheric temperatures in the hotspot are in the sample runsand therange of parameters for which equilib- range T 5 10 keV. The Fermi temperature falls below rium solutions can be found. ∼ − 10 keV only for a thin layer ( 0.01Z (ψ), from eq. 7) 0 ∼ at the top and a major fraction of the mound has Fermi 3.1 Modelling the magnetically supported temperature larger than 100keV. Comparison to the tem- accretion mound peratureprofiles obtained including heat transfer effects by Brown & Bildsten (1998), shows that the temperature in- We will assume a hydrogen poor plasma (µ = 2) being e side the mound remains lower than Fermi temperature at confined in the mound (Brown & Bildsten 1998). We re- greater depths. Thus modelling the confined accreted mat- strict the analysis to the gaseous state before ions form a teras a degenerate mound is appropriate. liquidphase.Theelectrostaticcouplingparameter(Γ)gives We have solved the G-S equation for different forms a rough estimate whether matter is solid (Γ >> 1), liquid of the mound height profile (e.g. Litwin et al. (2001); (Γ 1) or gas (Γ<1) (e.g. Litwin et al. (2001)) ≃ Hameury et al. (1983)) Γ = Z2e2(4πn)1/3 Z (ψ) = Z (1 ( ψ )2) (11) kBT 3 0 c − ψp Z2 ρ 1/3 108K ψ 1.1 (9) Z (ψ) = Z exp( 2 ) (12) ≃ (cid:18)A1/3(cid:19)(cid:18)108g cm−3(cid:19) (cid:18) T (cid:19) 0 c − ψp Hence mounds of base densities < 108 g cm−3 are con- Z (ψ) = Z (1 ( ψ )4) (13) sidered for this work. To determine the appropriate form 0 c − ψp of equation of state, we first check if the plasma is non- Themoundheightprofiledependsonthemassloadingfunc- relativistic (γ = 5) orrelativistic (γ = 4) byevaluating the tionattheaccretiondiscandredistributionofmatterinthe 3 3 adiabatic index (γ = dlnp) from theexpression of fermionic shock region for which at present there is no clear knowl- dlnρ pressure for degenerate electron gas (Chandrasekhar 1967). edge. We resort to evaluating the density by specifying the (cid:13)c 2011RAS,MNRAS000,??–?? 4 Mukherjee and Bhattacharya Mean ψ for mounds with Z above threshold (Z ) Density mound on polar cap 0.6 c max 1.0 Zc=75m Zc=76m 0.5 0.8 m n K 0.6 0.4 Height (z) i 0.4 ψnction()0.3 0.2 ux fu00..26 Zc=77m Zc=80m 0.0 n fl 0.0 0.2 0.4 0.6 0.8 1.0 Mea0.5 Radius in Km 0.4 Figure2. Theshapeoftheaccretionmoundplottedalongrand zaxisinequal scaletoshowtherealaspectratio.Themoundis 0.3 likeathinflatlayeronthepole. 0.2 0 50 100 1050 50 100 150 Iteration steps Field line plot for Z = 55m c 0.06 Figure 4. Mean ψ as a function of iteration steps for different moundheightsabovethestabilitythreshold.Themeanψisseen tooscillatebetweenmultiplestates.BeyondacertainZcitpasses m K 0.04 throughdifferentstates randomly. n height) i 0.02 mmaaxgnimetuicmfibealdseBde=ns1i0t1ie2sGle.sTshtheavnal1u0e8Zgccwma−s3ch(oassednistcoukseseedp Z ( in Sec. 3.1). For case (a) the total mass of the mound is 9 0.00 10−13M⊙ andmaximumbasedensityis 4.7 107 g∼cm−×3. 0.0 0.2 0.4 0.6 0.8 1.0 For case (b) the total mass of the m∼ound×is 2.13 Radius in Km 10−12M⊙ andmaximumbasedensityis 6.8 10∼7 g cm−×3. ∼ × Field line plot for Z = 70m Themoundisintheshapeofaflatthinlayeronthesurface c 0.08 ofthestar,confinedwithinthepolarcap(Fig. 2).Contours of ψ from the solution, which represent the magnetic field m 0.06 lines(asB ∇ψ=0) areplottedinFig. 3.Fromthefigure n K we see that·the field lines are bent to support the pressure ht) i 0.04 of the confined matter. The distortion is more in case (b). g Field lines are pushed outwards from the initial configura- ei h tion resulting in bunchingof field lines and increase in field Z ( 0.02 strength. 0.00 0.0 0.2 0.4 0.6 0.8 1.0 3.3 Valid parameter space for existence of Radius in Km solution Hameury et al. (1983) mention in their work that for the Figure 3. Solution for Zc = 55 m (top fig) and Zc = 70m configuration ofmoundheightprofilestheyhadconsidered, (bottom fig), cases (a) and (b) in the text. Solid lines are field no solution was found for field lower than a critical value. lines fromG-Ssolution. The dash-dotted linerepresents the top ofthemound. We observe the same for different values of magnetic field and different mound height profiles. For a fixed magnetic fieldwefindthatfortheparabolic profile(eq. 11),astable moundheight asasimplefunctionofψ,subjecttothecon- solution exists for only up to a maximum threshold mound straint : ρ 0 as r R , so that the mound is confined p height (Z ). For mounds higher than this, the ψ func- → → max within the polar cap. In most of the analysis we have used tion keeps oscillating between multiple states with closed eq. (11)whichhasrelativelyshallow gradientsandhelpsin magnetic loops during the iteration process and conver- speeding up thenumerical convergence. gence to an unique solution is not reached. For magnetic field 1012 G the maximum height of a mound for a sta- ∼ ble solution was found to be Z 70m. For a mound 3.2 Solutions from the GS-solver max ∼ higher than this threshold, Z = 75m, the ψ at intermedi- c Wediscuss here results from two sample runs ate steps pass through states similar in nature to states 1 Case (a) : Z =55 m and 2 as shown in Fig. 5. The mean ψ is seen to oscillate c Case (b) : Z =70 m between two states as in Fig. 4. For higher mounds the c withthemoundheightprofilespecifiedbyeq. (11)and branchesofmeanψ bifurcatetomultiplestatessimilar toa (cid:13)c 2011RAS,MNRAS000,??–?? Phase dependent view of Cyclotron lines 5 Field lines for state 1 Field lines for state 2 a length scale Z . Hence from eq. (7), eq. (10) and eq. max 0.15 (11) we get m ρ AZ3/2 ; p Z5/2 K0.10 ∼ max ∝ max eight in 0.05 ∇p≃B·∇B → Rpp ≃ ZBm2ax (15) H Z B4/7 max 0.00 Litwin et al. (20∝01) have shown that ballooninginstability 0.0 0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0 will disrupt the equilibria if ∆β > 7.8R /[(γ 1)Z (ψ)], Radius in Km p − 0 whereβ istheratioofplasmapressuretomagneticpressure Figure 5. The ψ function at two intermediate iteration steps [p/(B2/8π)].Usingp kadρ5/3(eq. 10)andρ AZm3/a2x(eq. ∼ ∼ (79th and 80th) of the GS-solver for a mound of height 75m. 7)wecanwritethestabilitycriterionobtainedbyLitwinas Closed magnetic loops are seen to form which indicate loss of 4 equilibria.Atdifferentiterationsteps,theψpassesthroughstates log (Z )> 5.1+ log B (16) 10 max − 7 10 similarto states 1 and 2 depicted here without reaching conver- gence. whichisveryclose totheobserveddependenceofZmax and B asobtainedfromournumericalsolutions.Thuslimitrep- resented by eq. (14) may result from ballooning typepres- Maximum mound height vs field sure driven instabilities where curvature of magnetic field 100.0 cannolongersupporttheplasmapressure(eq. 15,eq. 16) and theequilibrium solution cannot beobtained. Hencefor m 10.0 our analysis of the cyclotron line features in the following n section, we restrict ourselves to mounds of height less than imax 70m for a dipole field of 1012G. A detailed study of the Z 1.0 stability analysis of our solutions will be reported in a fu- turepublication (Mukherjee,Bhattacharyaand Mignonein preparation). 0.1 7 8 9 10 11 12 13 log (B) 10 4 CYCLOTRON RESONANCE SCATTERING Figure 6. Maximum height of the mound (Zmax) that can be FEATURES (CRSF) supportedbydifferentsurfacefieldstrengths.Moundheightpro- The major part of the mound forms an optically thick fileofeq. (11)isassumed.Thedashedlineshowsapower-lawfit mediumwithcyclotronlineformationtakingplaceinathin tothepoints.Validsolutionscanbeobtainedonlyforparameters belowtheline,abovewhichG-Scodedoes notconverge. layerlocated atthetopsurface.Thedepthofthelineform- ing region may be estimated as l Z z where Z (the 0 0 ∼ − mound height profile) is the top height of the mound at a pitch-forkdiagram.Formationofclosedmagneticloops,also given r. From the definition of optical depth and using eq. reportedinHameury et al.(1983),Payne& Melatos(2004) (7)forthedensitydistributionwefindtherelation between andPayne & Melatos(2007),appeartoindicatelossofequi- optical depth and thickness of the line forming region as libria. We have tried different simulations to check for the Aσ existenceofsolutionsbeyondthethreshold,e.g.higherreso- τ = l5/2 (17) µ m lutionrunsorimprovingtheinitialguesssolutionbystarting e p fromapreviouslyconvergedsolutionorincreasingtheradial where mp is proton mass and σ is the scattering cross sec- extent of thegrid. However stable solutions were not found tion. For Thomson scattering (σ = σTh), τ 1 occurs at ≃ for heights greater than thethreshold Z . a depth of 1.1 mm and for cyclotron resonance scatter- max ∼ Fig. 6 shows the maximum values of Zmax (eq. 11) ing (σ ∼ 105σTh) τ ≃1 occurs at 11.3 µm. Thus cyclotron for a given field up to which solutions exist. The maximum line formation takes place is a thin layer at the top of the allowed height (Z ) has a power-law dependenceon B mound. Variation in the local magnetic field at the top of max the mound is expected to cause variation in the cyclotron log (Z )= 3.676+0.461log (B) (14) 10 max − 10 spectra.Modellingofthecyclotronresonancescatteringfea- whereZ isinmetresandBinGauss.Asimilarmagnetic tures(hereafterCRSF)takingintoaccountthecontribution max fieldtoZ scalingwasobservedforsolutionwithdifferent of different parts of the mound to the line of sight is pre- max moundheightprofilese.gforeq. (12)wegetlog (Z )= sented in thefollowing section. 10 max 3.79+0.46log (B)andforeq. (13)wegetlog (Z )= − 10 10 max 3.61+0.45log (B) which indicates that this is a generic − 10 4.1 Modelling cyclotron spectra featureoftheGrad-Shafranovequationinthecurrentsetup. Such a scaling relation can be understood approxi- The emission profile from a point on the mound depends matelybycomparingthevariationofpressureandmagnetic on the strength and direction of the magnetic field and the field over different length scales which balance each other. angle between the emergent radiation and the local normal Lateral variation in pressure over scale R is balanced by to themound surface, which vary with position dueto cur- p tension from curvature in magnetic field which occurs over vature of the field lines and the shape of mound surface (cid:13)c 2011RAS,MNRAS000,??–?? 6 Mukherjee and Bhattacharya respectively.Wedivide themoundinto small area elements Light curve for single mound (Z= 45m) c (∆A in a (r,θ) grid in cylindrical coordinates on the 1.1 ri,θj polar cap) and find the local magnetic field vector and lo- cal normal to the mound for each element, assuming them 1.0 to constant over the area element. The resultant cyclotron ux spectrum is constructedbyintegrating theweighted contri- d fl 0.9 blinuetioonf soigfhetm(ihsesiroenaftferormlosa)ll parts of the mound towards a malise 0.8 01..9090 01..9090 r 0.98 0.98 o F = I(θ )∆A cosθ (18) N 0.97 0.97 αl ri,θj αl 0.7 0.96 0.96 i,j X 0.95 0.95 0 5 10 15 20 0 5 10 15 20 I is the normalised intensity at the mound surface and the 0.6 angleθ (eq. B6)istheanglebetweenthedirectionofemis- 0 60 120 180 240 300 360 αl Phase in degree sion at the mound surface and the local normal to the sur- face. For integrating the intensity over the mound we have chosen aradial grid with resolution equal toor higher than Figure 7. The light curve of one pole mound for Zc = 45 m, theresolution in theradial direction of the G-S solution. ηp =10◦ andi=30◦ (Appendix B).Insets showthespectra at the two rotation phases marked. No significant variation of line The cyclotron line energy is given by energy is seen. A detector resolution of 10% was assumed (eq. 20). n E E =nE √1 u 1 c0 sin2θ (19) n c0 − − 2 511keV αb (cid:18) (cid:18) (cid:19) (cid:19) where n = 1,2,3... is the order of the harmonic, Ec0 = 4.2 Results : cyclotron spectra from a single 11.6B12 in keV, θαb is the angle between the direction of mound emission andlocal magneticfield(eq. B5)andu=r /r,r s s beingtheSchwarzschildradius.Thefactor√1 ugivesthe Wehavecarried outsimulation runsfor thecyclotron spec- − tra using the solutions of the magnetic field obtained from gravitational redshift of theline. theGS-solver.Formostofouranalysisinthissectionweuse Eq. (19) is correct to second order in thesmall param- the solution with profile eq. (11). Although the shape and eter E /(511keV), and is adequate for the field strengths c0 natureoftheCRSFwillchangefordifferentprofilefunctions we consider in this work. In our studies we consider only wecandrawsomegeneralconclusionsaboutthedependence the effect of the fundamental n = 1 line. The accurate de- of thecyclotron spectra on the local magnetic field. pendenceofthewidthand depthof theCRSFontheangle We first consider the case of emission from a single θ can be obtained by solving radiative transfer equations αb hotspot at one of the poles. Fig. 7 shows the light curve inthemound.Thisisbeyondthescopeofthepresentwork from a mound of height Z = 45 m at a pole with incli- andwillbeaddressedinafuturepublication(Kumaretalin c nation η = 10◦ and los i = 30◦ (See Appendix B for preparation). For our present work we model the cyclotron p definitions of i and η ). The inset plots show the cyclotron featurebyaGaussianprofilewithlinecentresfromeq. (19) p spectra convolved with a Gaussian with f = 0.1 (eq. 20), andmodelthedependenceontheangleθ ofthelinewidth αb at two rotation phases 32◦ and 200◦. Although the magni- ∆E/E andtherelativedepthbyinterpolatingfromthere- c tude of the field at the top of the mound varies by 27% sults presented in Sch¨onherret al. (2007) for the slab 1-0 ∼ in this case (Fig. 8), the line centre of the CRSF shows geometry. less than 0.2 % change about a mean 9.6 keV. As the We incorporate effects of gravitational light bending ∼ continuumemissionisassumedtobeisotropicanduniform, (eq. B1) following the approximate formula given by and also since gravitational bendingredirects thelight rays Beloborodov (2002). For the intrinsic intensity profile we in directions well away from straight trajectories, all parts use a form I(θ )=I +I cosθ , but set I =0 for most αl 0 1 αl 1 of the mound contribute towards a given line of sight at all of our analysis. To simulate the finite energy resolution of phases. This gives a resultant averaged spectrum with very the detectors, we convolve the spectra with a normalised little phase dependenceof the CRSFfrom a single pole. Gaussian with the standard deviation a fraction f of the For mounds with large distortion of surface magnetic local energy field (e.g. for Z =70 m,thefield at theedgerises to 4.6 c ∼ 1 (E E′)2 times the original dipole value (see Fig. 8). We find two W(E,E′)= √2πσ exp − 2−σ2 ; σ∼fE′ (20) distinct CRSF fundamentals at two distinct energies (Fig. (cid:18) (cid:19) 9).Thefeatureat alowerenergy,forafield 1012 G,orig- ∼ Detectors currently used for observations in X-ray astron- inates near the centre of the mound while that at a higher omy usually have a energy resolution of 10%-20% (f energy arises in the high field regions near the periphery ≃ 0.1 0.2).Wecarryouttheaboveanalysisatdifferentphases of the mound (Fig. 8). We emphasize that this multiple- − ofrotationoftheneutronstartoperformaphasedependent featured spectrum is a result of the variation of the local studyofthespectra.Thenatureofthespectraalsodepends fieldstrengthanddoesnotrepresentmultipleharmonics,as on the relative orientation of the mound (inclination angle onlyn=1featureshavebeenincludedinourcomputation. η ) and the los (angle i) with respect to the spin axis of Theenergyratioofthesefeaturesmaythereforebearbitrary p the neutron star, which are treated as free parameters (see and not follow a harmonic relation. Appendix B for details on the geometry). Thelinecentresoftheindividualpeaksshowlittlevari- (cid:13)c 2011RAS,MNRAS000,??–?? Phase dependent view of Cyclotron lines 7 Variation of B at the top of the mound Light curve for single mound (Z= 60m) c 10 1.1 9 Zc=45m, Boffset= 8 1.0 x 8 Zc=50m, Boffset= 7 flu d 0.9 e 7 Zc=55m, Boffset= 6 alis 1.000 1.000 G) m 0.8 0.995 0.995 1210 6 Nor 0.990 0.990 / eld (B 5 Zc=60m, Boffset= 4 0.7 00..998805 00..998805 c fi 0.6 0 10 20 30 0 10 20 30 neti 4 0 60 120 180 240 300 360 ag Phase in degree M 3 Figure10. ThelightcurveofonepolemoundwithZc∼60m, 2 ηp = 10◦ and i = 30◦ (Appendix B). Insets show the spec- Zc=70m, Boffset= 0 tra at two the rotation phases marked. Multiple CRSF are seen withvariationoftherelativedepthwithphase,butnophasede- 1 pendence of the line energies. A detector resolution of 10% was assumed(eq. 20). 0 0.0 0.2 0.4 0.6 0.8 1.0 Radius in Km ationwithphasebuttherelativedepthofthepeaksdepend Figure 8. The variation of the strength of the magnetic field ontheviewinggeometryandthephaseangle.Fig. 10shows withradialdistanceatthetopofthemound,fordifferentZc.The the light curve and spectra (with f = 0.1 in eq. 20) from plotsareoffsetbyanamount(Boffset)forclarity.Themaximum a mound Z 60 m with η = 10◦ and i = 35◦, where magnetic fieldat the top isseveral times the surfacedipolefield c ∼ p therelativedepthoftheCRSFvarywithphase.Thedegree duetodistortionfrompressureofaccretedmatter. of their variation dependson the relative orientation of the pole and the los, occurring for a certain range of viewing Spectra from different mounds angles (25◦ i 40◦ for a pole at η =10◦ in the present p 1.10 ≤ ≤ case).TwodistinctCRSFareobservedforallmoundsofap- preciablemagneticfielddistortionatthetopandallviewing Zc=50m geometries, which shows that it is a generic feature of the 1.08 cyclotron lines originating from themound. Thecyclotronfeaturesaredominatedbythefieldstruc- turesneartheperipheryofthemoundastheseregionshave 1.06 alargeremittingarea.ForG-Ssolutionswithmoundheight profile as in eq. (11), the regions of higher magnetic field x u Z=55m d fl c distortion occur near the edge of the polar cap. Hence we e observetwoCRSFofcomparabledepthsforallmoundswith alis 1.04 profile eq. (11). However, if we use an exponential profile m or (eq. 12)whichfallsoffsharplywithr,regionsoflargerfield N distortion occur much closer to the axis (Fig. 11), and the 1.02 Zc=60m CRSFisdominatedbytheundistortedfieldintheouterre- gion. A shallow feature at 18 keV is contributed by the ∼ higher field regions nearer to the axis. This shows that the 1.00 CRSFisheavilyinfluencedbythegeometryanddistribution of thelocal field as well as thestructureof themound. Z=70m c 0.98 0 20 40 60 4.3 Effects of finite energy resolution of detectors Energy in KeV The finite energy resolution of the detector can make the Figure9. CRSFfrommoundsofdifferentheights.Thespectrum twoabsorptionfeaturesindistinguishableifcloselyplaced.In for each mound shows two CRSF fundamentals at two distinct Fig. 12weshowtheeffectofGaussianconvolutionwithdif- energies, corresponding to the undistorted dipolar field and the ferentvaluesoff (eq. 20).Forf =0.2,thetwopeaksarein- regionwithlargedistortionofthemagneticfieldrespectively(see distinguishable.Formoundsoflowerheight(e.g.Z =45m) c Fig. 8). the two CRSF become indistinguishable at a much smaller f.Thusthefiniteenergyresolutionofthedetectorcanoften mask the effects of the internal magnetic field structure on the CRSF. However for mounds of higher height, the two (cid:13)c 2011RAS,MNRAS000,??–?? 8 Mukherjee and Bhattacharya Exponential profile with Z = 55m Two antipodal mounds (45m & 52m) c 0.06 1.01 m 1.00 K 0.04 n x ht) i flu 0.99 Z (heig 0.02 alised 0.98 1.00 1.00 m 0.99 0.99 0.00 Nor 0.97 Ec=9.4KeV Ec=11.4KeV 0.0 0.2 0.4 0.6 0.8 1.0 0.98 0.98 Radius in Km 0.96 0.97 0.97 Z=55m, Exponential profile 0 10 20 30 0 10 20 30 1.00 c 0.95 0 60 120 180 240 300 360 0.99 Phase in degree x u malised Fl 00..9978 4F5igmuraend1532. mT,hweitlhigηhpt=cu3rv0e◦faonrdt1w5o0◦anretis-ppeocdtaivlemly,oaunnddsi,=Zc30∼◦ or (Appendix B).Insetsshowthespectraatthetworotationphases N 0.96 marked.TheCRSFshows20%changeinlineenergywithphase. Adetector resolutionof15%wasassumed(eq. 20). 0.95 0 10 20 30 Energy in KeV CRSF have centres far enough to be distinguishable even Figure 11. Magnetic field line structure (upper panel) for a with f =0.2. moundwithZc∼55mwithanexponentialmoundheightprofile eq. (12)andthespectrum emitted(lowerpanel). Thespectrum is dominated by area elements near the polar cap edge (corre- sponding to the region of undistorted field in this case) as they 4.4 Emission from two anti-podal poles span larger areas. The shallow feature at ∼18 keV corresponds totheregionoflargefielddistortionclosertotheaxis. Many systems show variation of the cyclotron energy with the phase of rotation of the neutron star e.g. ∆E /E 10%forVelaX-1(Kreykenbohmet al.2002), cyc cyc Effect of Gaussian convolution ∼ 20% for 4U0115+63 (Heindlet al. 2004; Baushev 2009), 1.15 ∼ 25% for Her X-1 (Klochkov et al. 2008b), GX 301- ∼ f=0.2 2 (Kreykenbohmet al. 2004), and 30% for Cen X-3 ∼ (Suchyet al. 2008; Burderi et al. 2000). From our simula- tions of CRSF for emission from a single hotspot we find 1.10 f=0.15 very little variation of the line energy with rotation phase. However if emission is considered from two antipodal host- potsatoppositepoleswithslightdifferenceinmoundheight dueto asymmetric accretion, then the line centre of there- ed flux f=0.1 ssuhlotwanstthCeRliSgFhtschuorwveafosrtraonngeeurtrpohnassteardewpiethndtewnocea.nFtii-gp.od1a3l malis 1.05 poles at ηp ∼30◦ and 150◦, having mounds of height Zc ∼ Nor f=0.05 45 m and 52 m respectively, and an observer at inclination i = 80◦. The simulation is carried out by assuming uni- form and isotropic continuum intensity normalised to 1 at both poles, which may not be valid in reality due to differ- 1.00 No convolution ence in accretion rates at the two poles. However for small differencesin themoundheights considered here, weignore thedifferentialluminosityeffectanddrawsomegeneralcon- clusions about the behaviour of the cyclotron line energies. The spectra are convolved with a Gaussian function with 0.95 0 5 10 15 20 25 30 f =0.15 (eq. 20). Energy in KeV Thelightcurveisnotsinusoidalunlikethecaseofasin- glepole(Fig. 7andFig. 10).ThelineenergyoftheCRSF Figure 12. The effect of finite detector resolutionmodelled by varies by ∆E 20% over a full rotation cycle. The spec- c convolvingthespectrawithaGaussianfunction.Resultsfordif- ∼ trum from the pole nearer to the los dominates the CRSF. ferentvaluesoff (eq. 20)areshown.ThedoubleCRSFnatureof Theobserved variation inthelineenergyof coursedepends thespectradisappearforf ≥0.2correspondingtoa20%energy on the viewing geometry defined by the angles i and η . A resolutionofdetector. p proper evaluation of the spectra for a two pole case would require the knowledge of the accretion rate and local tem- perature of the hotspots, but it may be concluded that in (cid:13)c 2011RAS,MNRAS000,??–?? Phase dependent view of Cyclotron lines 9 Spectra for Z = 60m, η = 100, i=700 (ii) Multiple, anharmonic CRSF : From our simula- c p 1.000 tionsweobserveMultipleCSRFfundamentalsfromasingle x 0.995 mound with considerable field distortion (see Fig. 9, Fig. ed flu 0.990 1m0o)n.iAcsrweahlicshpewctilrluamddwitlol atlhseoccoomntpalienximtyulotfiptlheeCsRpSecFtrhuamr-. s mali 0.985 SeveralHMXBsshowmultipleCRSFwherethelinecentres or have anharmonic separation e.g. Ec 22 keV, 47 keV for N 0.980 4U 1538-52 (Rodes-Roca et al. 2009),∼ 23keV,51keVfor ∼ 0.975 VelaX-1(Kreykenbohmet al.2002), 45keV,100keVfor ∼ 0 10 20 30 A 0535+26(Kendziorra et al. 1994; Caballero et al. 2007), Energy in KeV 26 keV,49 keV,74keV for V 0332+53 (Makishima et al. ∼ 1990; Pottschmidt et al. 2005), 36 keV, 63 keV for EXO Figure14. TheCRSFfromamoundofheightZc=60mwith 2030+375 (Reig & Coe 1999; Kl∼ochkov et al. 2008a). More mound height profile eq. 11, ηp = 10◦ and i = 70◦ (see Ap- on the properties of these sources can be found in the re- pendix B). The undistorted field produces the small feature at views by Mihara et al. (2007) and Lutovinov & Tsygankov ∼8keV.Thedistortedfieldnearthepolarcapedgeisthedomi- (2008). nantcontributortotheCRSF,producingadipat∼18keV.This Anharmoniclinespacingduetointrinsicnon-dipolarfield shows that field inferredfrom CRSFline energy may not depict the surface dipole field, but the field distorted due to accreted has been discussed by Nishimura (2005). Our results show matter. that strong non-dipolar structure can be generated in the accretion process itself. (iii) Effect of detector resolution : The detectability thepresenceofmultiplehotspots,theCRSFlineenergywill of multiple cyclotron features depend on the energy resolu- show a significant phase dependence. tion of the detector. For detectors with poor energy resolu- tion, the multiple fundamental features will be masked as shown in Fig. (12). However for mounds with large field distortion(Z 60 70m)multipleabsorptionfeatureswill c ∼ − 5 DISCUSSION AND CONCLUSIONS still be observed with detectors that haveenergy resolution of 20%. We have modelled the structure of the accretion mound ≤ (iv) CRSF phase dependence - one pole : We find by solving for the static magnetic equilibria described by no appreciable spin phase dependence of the line energy of theGrad-Shafranovequation. Formoundsof total accreted mass 10−12M⊙ (Sec. 3.2) there is appreciable distortion a given CRSF from a single mound (Fig. 7) despite there ∼ being considerable local variation of magnetic field. This ofthemagneticfield.Wehavefoundthatforagivensurface may be attributed to the fact that all parts of the mound field strength, stable solutions to the G-S equation are not contribute towards a given line of sight at any phase due found for mounds of height higher than a threshold (Sec. to effects of gravitational light bending and assumption of 3.3). Beyond thethreshold,field lines with closed loop con- isotropic local emission. However we find that the relative figurations are formed, and may indicate the onset of pres- depth of themultiple CRSF dependson phase and viewing suredrivenMHDinstabilities.Sowerestrictouranalysisto geometry as shown in Fig. 10. mound heights for which equilibrium solution is obtained. (v) CRSFphasedependence-twopole:Lineenergy Using the mound structure obtained from the Grad- variation is however observed if emission from both poles Shafranov equation, we have simulated the cyclotron spec- with slightly different mound heights are considered (Sec. tra which will originate from its mound surface. We have 4.4). For emission from two anti-podal mounds of height integratedtheemissionfromdifferentpartsofthemoundto Z 45 m and 52 m the CRSF line energy varies by find the resultant spectra and have performed a phase de- c ∼ ∼ 20% during one spin cycle, similar to what is observed for pendentstudyof thespectra. Wehaveassumed aGaussian manysources.Thusweconcludethatemissionfrommultiple modelfortheCRSFfundamentalandhaveincorporatedef- accretion mounds will result in strong variation of the line fects of gravitational bending of light and finite energy res- energy of theCRSFwith spin phase. olutionoftheX-raydetectors(Sec. 4.1).Fromouranalysis (vi) Dependence on mound structure : CRSF are ofthephasedependentspectra(Sec. 4.2,Sec. 4.3andSec. dominated by field structure at the mound periphery due 4.4) we can draw some general conclusions : tothelargerphysicalareaofthisregion.Differentstructure (i) CRSF need not represent the dipole field :The anddensitydistributionofthemoundwouldcausedifferent CRSF are heavily influenced by the distortion of the local distributionsofthelocal field.Fig. 9andFig. 11showthe field caused by the pressure of the confined accreted mat- differenceintheCRSFfromamoundofthesamemaximum ter. Theline energy of theCRSFmay not always represent height(Z =55m)butdifferentmoundheightprofiles.The c the intrinsic dipole magnetic field. As shown in Fig. 14, strong dependence of the spectra on the structure and size thesmall dip at 8keV is theredshifted cyclotron feature of the mound suggests that with variation in the accretion due to the surfa∼ce field of 1012 G, which may not be ob- rate, the observed line profiles and energies may change, servedinpracticeinthepresenceofnoiseandpoordetector contributingto a luminosity dependenceof thespectrum. resolution. The CRSF at 18 keV resulting from the dis- Some sources e.g. Her X-1 (Staubert et al. 2007) and A ∼ torted field will be the dominant feature in the spectrum. 0535+26(Klochkov et al.2011),showapositivecorrelation For mounds of height Z 70 m the maximum distortion betweenluminosityandlineenergy.Inourpicture,thiscan c ∼ in thefield can be as large as 4.6 times thedipole value. be attributed to the presence of a strong non-dipolar com- (cid:13)c 2011RAS,MNRAS000,??–?? 10 Mukherjee and Bhattacharya ponentinthefieldduetoanincrease inmassofthemound toAndreaMignoneforhishelpandsuggestionsinsettingup resultingfromanincreaseinaccretionrate.Oursimulations simulation runs with PLUTO. We thank Dr. Kandaswamy show that for a small change of mound size we have an ap- Subramanian,Dr.RanjeevMisraandSandeepKumarfrom preciable change in maximum magnetic field at the top of IUCAA, for useful discussions and suggestions during the themound(Fig. 8),e.g.between Z of55mand60m,the work, and IUCAA HPC team for their help in using the c maximum field at thetop of themoundchanges by 36%. IUCAA HPC where most of the numerical computations ∼ Since we work in the static limit and do not consider ac- were carried out. We also would like to thank Dr Ru¨diger cretion rate as a parameter in our simulations, we cannot Staubertforhissuggestions in improvingthetext,andalso directly probe the luminosity dependence of the spectrum. theanonymousrefereeforhisvaluablesuggestionsandcom- However we can conclude that small changes in height of ments. themoundcanresult insignificantchangesinthemagnetic field inferred from theCRSF. (vii) Effect of anisotropic continuum : Wehavealso APPENDIX A: SOLVING THE carried out runs with mildly anisotropic continuum inten- GRAD-SHAFRANOV EQUATION sity profiles : I(θ ) = I (1+cosθ ) and have found no αl 0 αl A1 The numerical algorithm appreciable phase dependence of the line energy of a one pole CRSF although the percentage modulation of the flux We express the G-S equation in non-dimensional form by was larger. However if continuum is highly anisotropic, e.g. scaling the physical parameters with L =R =1km, ρ = 0 p 0 I = (sin2θαl/(2θαl))2 then the resultant spectra are found 106 g cm−3, B0 = 1012 G. The flux function is normalised to be phase dependent. However such strong beaming may to the value at r = R : u = ψ/ψ (ψ = 1B R2). We p p p 2 0 p not be realistic in thecontext of HMXBpulsars. performavariabletransformation x=r2 andsolvetheG-S (viii) Screening effect from overlying column : Ef- equation on a (x y) grid which are the normalised (r2 fects of partial screening of the mound by an atmosphere z) coordinates. Th−e result is then transformed back to th−e of accreted matter have not been considered in our simu- (r z)gridbysplineinterpolation.Withthistransformation lations. An optically thick blanket of settling plasma can an−d scaling, theG-S equation for an adiabatic gas reads as screen a fraction of the mound depending on the viewing ∂2u ∂2u dY (u) geometry. The effect of screening can be estimated approx- 4x + = x(Y y)α 0 (A1) imately by using the velocity profile of the settling flow ∂x2 ∂y2 −C 0− du nfreoumtroBnecsktaerr&ofWmoalsffs (2010.74)M:⊙va(zn)d =rad0iu.4s9c 1z0/kRmp,, mfoerana mwhoeurnedChe=ight16pπrAofiglLe10e+xαp/rBes02s,edαin=scaγl−e1d1 flaunxdcoYo0r(dui)naistetuh.e plasma temperature ∼10keV and magnetic∼pfield 1012G. Weadopt thefollowing set of boundary conditions : Foran accretion rate∼M˙ 1016gs−1 and amound∼of height ∼ (i) ψ=ψ atz=0 r.ψ = 1B r2 istheinitialguessψ ZTcho∼ms7o0nmop(tmicaaxlimdeupmthaalsloτwed1Z.5c7for1e0q−.5ℓ,11w)hewree gℓeits tthhee which is thedψ for a un∀iformdfield2B0=B0zˆ(approximation ∼ × ofdipolarfieldinthepolarcapregion).Thisistheline-tying lineelementalongthelosthroughtheaccretioncolumn.For ahotspotofradiusR 105cm,weseethatregionsnearthe condition of the field at thebase (Hameury et al. 1983). p ∼ (ii) ψ=ψ at r=0 z. axis will be optically thick along a radial line from the axis d ∀ (iii) ψ = ψ at r z. Field is undistorted for (in local cylindrical geometry of themound).The degree of d → ∞ ∀ r >> R . Ideally one should keep the r = r boundary screeningofdifferentpartsofthemoundwillvarywithspin p max at infinity to fully capture the distortion of the field lines phase, resulting in phase dependentspectra. resulting from the lateral pressure of the confined plasma. Thusweconcludefromthisworkthatthedistortionsin To implement this numerically is impractical. In our GS- local magnetic field due to confinement of accreted matter solver the boundary along the radial direction is chosen at at the polar cap has considerable influence on the phase r =r whereψ=ψ isfixedastheboundarycondition. max p d resolvedspectrafromaccretingbinaryneutronstarsystems. The mound height profiles considered for the GS-solver fall Current observations often have poor count statistics and off with radius vanishing at r =r . Thus at r =r bound- p p involvesaveraginginphaseduetowhichmanyofthedetails ary there is very little deviation from initial unloaded field inthespectramaybelost.FuturemissionslikeASTROSAT configuration. Possible limitations of this are discussed in (Agarwal et. al. 2005 , Koteswara Rao et al. (2009)) with Sec. 3.3. higher sensitivity and better energy resolution can help us (iv) ψ = ψ at z . Field is undistorted for z . d analyse the spectra in greater detail and provide us clues Setting boundary he→igh∞t (z ) too close to the max→imu∞m top to thenature, shape and geometry of the accretion mound. height of the mound (Z ) affects the solution and gives in- c A more detailed analysis including dynamic simulations of correctresult.Itisseenthatsettingz 1.5Z issufficient top c the emitting region and full radiative transfer would help for the solution to be stable. Further c≥hange in boundary usinvestigateeffectslikelineenergy-luminositycorrelation, height (z ) does not change thesolution significantly. top anharmonic separation of the CRSF energies, asymmetric Weadoptanumericalschemesimilartotheonefollowed shape of the CRSFetc. by Mouschovias (1974) and Payne& Melatos (2004). The density is determined from the mound height profile and eq. (7). A two state iterative scheme using an outer under 6 ACKNOWLEDGEMENT relaxation scheme and an innerSuccessive OverRelaxation We thank CSIR India for Junior Research fellow grant, (SOR) method using Chebyshev acceleration (Press et al. award no 09/545(0034)/2009-EMR-I. We are very thankful 1993)isadoptedtotacklethenon-linearsourcetermonthe (cid:13)c 2011RAS,MNRAS000,??–??