Thermal conductivity and phase separation of the crust of accreting neutron stars
aa r X i v : . [ a s t r o - ph ] J a n Thermal conductivity and phase separation of the crust of accreting neutron stars
C. J. Horowitz ∗ and O. L. Caballero † Department of Physics and Nuclear Theory Center, Indiana University, Bloomington, IN 47405
D. K. Berry ‡ University Information Technology Services, Indiana University, Bloomington, IN 47408 (Dated: December 1, 2018)Recently, crust cooling times have been measured for neutron stars after extended outbursts.These observations are very sensitive to the thermal conductivity κ of the crust and strongly suggestthat κ is large. We perform molecular dynamics simulations of the structure of the crust of anaccreting neutron star using a complex composition that includes many impurities. The compositioncomes from simulations of rapid proton capture nucleosynthesys followed by electron captures. Wefind that the thermal conductivity is reduced by impurity scattering. In addition, we find phaseseparation. Some impurities with low atomic number Z are concentrated in a subregion of thesimulation volume. For our composition, the solid crust must separate into regions of differentcompositions. This could lead to an asymmetric star with a quadrupole deformation. Observationsof crust cooling can constrain impurity concentrations. PACS numbers: 97.60.Jd, 26.60.+c, 97.80.Jp, 26.50.+x
I. INTRODUCTION
What is the thermal conductivity of the crust of a neu-tron star? Recently the cooling of two neutron stars hasbeen observed after extended outbursts [1, 2]. These out-bursts heat the stars’ crusts out of equilibrium and thenthe cooling time is measured as the crusts return to equi-librium. The surface temperature of the neutron star inKS 1731-260 decreased with an exponential time scaleof 325 ±
100 days while MXB 1659-29 has a time scaleof 505 ±
59 days [2]. Comparing these observations, ofrapid cooling, to calculations by Rutledge et al. [3] andShternin et al. [4] strongly suggest that the crust has ahigh thermal conductivity. This would be expected if thecrust is a regular crystal.In contrast, a low crust thermal conductivity, thatwould be expected if the crust is an amorphous solid,could help explain superburst ignition. Superbursts arevery energetic X-ray bursts from accreting neutron starsthat are thought to involve the unstable thermonuclearburning of carbon [5, 6]. However, some simulations donot reproduce the conditions needed for carbon ignitionbecause they have too low temperatures [7]. A low ther-mal conductivity could better insulate the outer crustand allow higher carbon ignition temperatures.The thermal conductivity is dominated by heat con-duction by electrons and this is limited by electron-ionscattering [8]. Therefore in this paper, we present molec-ular dynamics simulations of the crust in order to calcu-late electron-ion scattering. We include many impuritiesbased on results of a rapid proton capture nucleosynthesis ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] simulation [9] followed by calculations of electron capture[10]. We find a regular crystal structure. We do not findan amorphous phase. We calculate the static structurefactor S ( q ), that describes electron-ion scattering, andfrom S ( q ) we determine the thermal conductivity. Impu-rities can limit the thermal conductivity. If the impuri-ties are weakly correlated than their effect on the thermalconductivity can be described by an impurity parameter Q [11], Q = (∆ Z ) = h Z i − h Z i . (1)This depends on the dispersion in the charge Z of eachion. The rp process ash composition of ref. [10] and ref.[12] has a relatively large value of Q = 38 .
9. Impurityscattering can be important at low temperatures wherethere is small scattering from thermal fluctuations. Notethat ref. [11] assumes the impurities are weakly corre-lated. If there are important correlations among the im-purities, for example if there is a tendency for low Z ionsto cluster together instead of being distributed at randomthroughout the lattice, then the effects of impurities onthe thermal conductivity could be different from what iscalculated in ref. [11]. In this paper we perform MD sim-ulations to study the distribution of impurities and theireffect on the conductivity.If the thermal conductivity is high, one may need ad-ditional heat sources in the crust in order to explain su-perburst ignition. Although Gupta et al. [10] find someheating from electron captures to excited nuclear states,simple nuclear structure properties may provide a natu-ral limit to the total heating from electron captures [13].Horowitz et al. [14] find additional heating from fusionof neutron rich light nuclei such as O+ O at densitiesnear 10 g/cm . These fusion reactions are an importantarea for future work. Alternatively chemical separationwith freezing, that was found in ref. [12], could enrichthe neutron star ocean with low Z elements and make it TABLE I: Abundance y z (by number) of chemical element Z . Z Abundance8 0.030110 0.011612 0.002314 0.002315 0.002320 0.004622 0.081024 0.071826 0.101927 0.002328 0.076430 0.085632 0.011633 0.125034 0.386636 0.002347 0.0023 easier for superburst ignitition.In section II we describe our molecular dynamics sim-ulations and the calculation of the thermal conductivity.Results for the structure of the crust, the static structurefactor S ( q ), and the thermal conductivity are presentedin section III. We conclude in section IV. II. MOLECULAR DYNAMICS SIMULATIONS
In this section we describe our molecular dynamics sim-ulations and how we calculate the thermal conductivity.We begin with a discussion of our initial composition.
A. Compositon
Our model for the composition of the crust is the sameas was used in previous work on chemical separationwhen the crust freezes [12]. Schatz et al. have calcu-lated the rapid proton capture (rp) process of hydrogenburning on the surface of an accreting neutron star [9],see also [15]. This produces a variety of nuclei up tomass A ≈ . × g/cm , has forty % ofthe ions with atomic number Z = 34, while an additional10% have Z = 33. The remaining 50% have a range oflower Z from 8 to 32. In particular about 3% is O and1% Ne. This Gupta et al. composition [10] is listed inTable I. In general, nuclei at this depth in the crust areexpected to be neutron rich because of electron capture.Material accretes into a liquid ocean. As the den-sity increases near the bottom of the ocean, the materialfreezes. However we found chemical separation when thecomplex rp ash mixture freezes [12]. The ocean is greatly enriched in low Z elements compared to the newly formedsolid. What does chemical separation mean for the struc-ture of the crust? Perhaps the most conservative possi-bility is the following steady state scenario. We assumematerial accretes at a constant rate. Initially, chemicalseparation enriches the ocean in low Z elements. Eventu-ally the ocean becomes so enriched that the compositionof low Z material in the newly forming solid is equalto that in the accreting material. The system reachesa steady state. The rate of low Z material accretinginto the ocean is equal to the rate freezing out (modulonuclear reactions). The sole effect of chemical separa-tion is to greatly enrich the ocean in low Z material.If we assume steady state, the composition of the crustwill be the same as that of the original accreting mate-rial. Therefore, in this paper we perform MD simulationsto determine the structure and thermal conductivity ofcrust with the original Gupta et al. rp ash composition. B. Simulations of Crust Structure
In order to calculate the thermal conductivity of a mul-ticomponent system one needs to understand its state.Monte Carlo simulations [16] of the freezing of a classicalone component plasma (OCP) indicate that it can freezeinto imperfect body centered cubic (bcc) or face-centeredcubic (fcc) microcrystals. Unfortunetly not much hasbeen published on the freezing of a multi-componentplasma (MCP). There are many possibilities for the stateof a cold MCP [17]. It can be a regular MCP lattice; ormicrocrystals; or an amorphous, uniformly mixed struc-ture; or a lattice of one phase with random admixture ofother ions; or even an ensemble of phase separated do-mains. We perform classical MD simulations to explorethe state of our MCP solid.The electrons form a very degenerate relativistic elec-tron gas that slightly screens the interaction betweenions. We assume the potential v ij ( r ) between the ithand jth ion is, v ij ( r ) = Z i Z j e r e − r/λ e , (2)where r is the distance between ions and the electronscreening length is λ e = π / / [2 e (3 π n e ) / ]. Here n e isthe electron density. Note that we do not expect our re-sults to be very sensitive to the electron screening length.For example, the OCP melting point that we found inref. [12], using a finite λ e , agrees well with the result for λ e = ∞ .To characterize our simulations , we define an averageCoulomb coupling parameter Γ for the MCP,Γ = h Z / ih Z i / e aT , (3)where the mean ion sphere radius is a = (3 / πn ) / and n = n e / h Z i is the ion density. The OCP freezes at Γ =175. In ref. [12] we found that the impurities in our MCPlowered the melting temperature until Γ = 247. Finally,we can measure time in our simulation in units of oneover an average plasma frequency ω p , ω p = (cid:16)X j Z j πe x j nM j (cid:17) / , (4)where M j is the average mass of ions with charge Z j and abundance x j (by number). Note that there willbe quantum corrections to our classical simulations fortemperatures significantly below the plasma frequency. C. Thermal conductivity
The thermal conductivity κ has been discussed byPotekhin et al. [8]. We assume κ is dominated by heatcarried by electrons [8], κ = π k B T n e m ∗ e ν , (5)where the effective electron mass is m ∗ e = ǫ F = ( k F + m e ) / with k F the electron Fermi momentum and m e the electron mass. The electron collision frequency ν isassumed to be dominated by electron-ion collisions [8], ν = 43 π h Z i ǫ F α Λ . (6)Here α is the fine structure constant and Λ is theCoulomb logarithm that describes electron-ion collisions[8], Λ = Z k F q dqqǫ ( q, S ′ ( q )(1 − q k F ) . (7)Here ǫ ( q,
0) is the dielectric function due to degeneraterelativistic electrons, [18], see Eq. 2.3 of [19]. Note thatfor simplicity we neglect second and higher Born correc-tions to electron ion scattering in Eq. 7, see for exam-ple [19]. We are interested in the difference in thermalconductivity for different solid structures. Second andhigher Born corrections should be the same for the dif-ferent structures. Finally, the lower limit q in Eq. 7 is0 in a liquid phase and q = (6 π n ) / in a crystal phase[8].The static structure factor S ( q ) describes electron-ion scattering. We calculate S ( q ) directly as a density-density correlation function using trajectories from ourMD simulations, S ( q ) = h ρ ∗ ( q ) ρ ( q ) i − |h ρ ( q ) i| . (8)Here the charge density ρ ( q ) is, ρ ( q ) = 1 √ N N X i =1 Z i h Z i e i q · r i , (9) with N the number of ions in the simulation and Z i , r i are the charge and location of the ith ion. We evaluatethe thermal average in Eq. 8 as a time average duringour MD simulations.The static structure factor S ( q ) includes both Braggscattering contributions from the whole crystal lattice S bragg ( q ) and inelastic excitation contributions S ′ ( q ) [8], S ( q ) = S ′ ( q ) + S bragg ( q ) . (10)The Bragg contribution is a series of delta functions atmomenta related to one over the lattice spacing. Thisdescribes Bragg scattering and helps determine the elec-tron band structure. It does not limit the electron meanfree path. Instead the mean free path and thermal con-ductivity are determined by S ′ ( q ).Our MD simulations are classical. Unfortunately thisclassical approximation makes the separation of S ( q ) into S ′ and S bragg somewhat ambiguous. We approximate S ′ ( q ) with a simple numerical filter applied to S ( q ). Thefilter removes delta function like contributions to S ( q )that have a very rapid q dependence, and also removesnumerical noise. This is discussed further in Section III. III. RESULTS
To explore possible states for the multicomponentplasma we perform two molecular dynamics simulations.The initial conditions of these simulations are similar tothose in [12]. The composition is indicated in Table I. Westart by freezing a very small system of 432 ions. Herethe ions were started with random initial conditions at ahigh temperature T and T was reduced in stages (by re-scaling velocities) until the system freezes. For the firstsimulation run, called rpcrust-01b in Table II, we placefour copies of this 432 ion solid in a larger simulationvolume along with four copies of a 432 ion liquid con-figuration. This 3456 ion configuration is evolved at alower temperature until the whole system freezes. Next,we evolve the 3456 ion solid at a reference density of n = 7 . × − fm − (or 1 × g/cm ) for a totalsimulation time of 2 . × fm/c (8 . × ω − p ). Thetemperature was started at 0.325 MeV and slowly de-creased to a small value by the end of this time. Thedensity and initial temperature correspond to Γ = 261 . t = 25 fm/c for a total of9 . × steps. This took about 2 months on a singlespecial purpose MDGRAPE-2 [21] board. Next, this lowtemperature configuration was reheated to T = 0 . . × fm/c. The total time was4 × fm/c. This somewhat complicated procedure wasdone for historical reasons. It does allow plenty of timefor ions to diffuse throughout the simulation volume.Note that at our artificially high reference density (10 g/cm ) free neutrons will be present. However, we areprimarily interested in lower densities with out free neu-trons. Our results can be scaled to other densities and TABLE II: Computer Simulations. The start time is t i , thefinish time t f , N is the number of ions, and the temperatureis T . Each simulation is at a density n = 7 . × − fm − (1 × g/cm ). Note that one over the plasma frequency is1 /ω p = 270 fm/c.Run N t i (fm/c) t f (fm/c) T (MeV)rpcrust-01b 3456 0 1 . × × × × × . × . × temperatures such that the Coulomb parameter Γ re-mains the same, see below. Furthermore, although wequote all simulation times in fm/c, the times can be ex-pressed in terms of one over the average plasma frequencyusing 1 /ω p = 270 fm/c.The initial configuration for run rpcrust-01b, see TableII, is shown in Fig. 1. The system is seen to be composedof two micro-crystals of different orientations. This issimilar to the micro-crystals found in ref. [16] upon freez-ing a one component plasma. In Fig. 1 we highlight thepositions of the O ions (as small red spheres). Theseions are located both in the crystal planes and in betweenthem. The O ions are not spread uniformly throughoutthe volume but there is a tendency for them to cluster.This will be discussed in more detail below.This configuration was then reheated to T = 0 . . × fm/c. The final config-uration of run rpcrust-01b is shown in Fig. 2. The twomicro-crystals of different orientation are now gone. Thesystem has managed to anneal into a single crystal witha single orientation. This suggests that micro-crystalscould be an artifact of computer simulations of limitedsize and duration. It also suggests that neutron star crustcould be formed with relatively large domain sizes.Figure 2 shows that O ions and other low Z impuri-ties are enhanced in regions on the left and right of thesimulation volume. Because of the periodic boundaryconditions this actually corresponds to a single region.We conclude that this complex mixture does not form asingle uniform solid phase. Instead it separates into twosolid phases. One phase is enriched in high Z ions andthe other phase is enriched in low Z ions.To study this further we have performed another sim-ulation labeled rpcrust-05 in Table II. The starting pointwas similar to run rpcrust-01b with eight copies of a 432ion configuration placed into a larger simulation volume.This 3456 ion configuration was evolved for 2 . × fm/c as the temperature was slowly decreased from 0 . T = 0 . T = 0 . T = 0 . T = 0 .
1, 0.2, and 0.3 MeV runs was 3 . × fm/c.The final configuration of run rpcrust-5 is shown in Fig.3. The system involves only a single body-centered cubic(bcc) crystal. However O and other low Z ions are notuniformly distributed. Instead they are strongly enrichedin a local region. This is indicated in Fig. 4 that showsthe radial distribution function g ( r ) for run rpcrust-05at a temperature T = 0 . g ( r ) for O-O is seen to be larger than one overa range of moderate distances r . This shows that the Oions are concentrated in a localized sub-volume. We con-clude that the complex rp ash mixture does not form asingle solid phase. Instead, for this composition, the neu-tron star crust must be composed of two or more regionsof different compositions. This disproves our steady stateassumption. There appears to be no composition of theliquid ocean, no matter how enriched in low Z ions, thatallows a uniform solid phase to form.These multiple regions of the crust with different com-positions may be very important for the structure of theneutron star. For example, if the phases are not dis-tributed uniformly, this could lead to a mass quadru-ple moment that might radiate gravitational waves [23].This nonuniform distribution of phases could arise froman anisotropic temperature because phase separation istemperature dependent.Finally for comparison we have also performed a onecomponent plasma simulation, see run OCP in Table II,where each ion has a charge Z = 29 . FIG. 1: (Color on line) Configuration of the 3456 ion mixturein run rpcrust-01b at the start of the simulation. The smallred spheres show the positions of O ions, while ions of aboveaverage Z are shown as larger blue spheres. Finally, ions ofbelow average Z (except for O) are shown as small whitespheres. The left and right halves of the figure show twomicro-crystals of different orientations. FIG. 2: (Color on line) Configuration of the 3456 ion mixturein run rpcrust-01b after a simulation time of 1 . × fm/c.The small red spheres show the positions of O ions, whileions of above average Z are shown as larger blue spheres.Finally, ions of below average Z (except for O) are shown assmall white spheres.FIG. 3: (Color on line) Configuration of the 3456 ion mixturein run rpcrust-05 after a simulation time of 1 . × fm/c.The small red spheres show the positions of O ions, whileions of above average Z are shown as larger blue spheres.Finally, ions of below average Z (except for O) are shownas small white spheres. The O concentration is seen to beenhanced in a sub-region to the right of center.
A. Static Structure Factor
We calculate the static structure factor S ( q ) from thedensity-density correlation function, Eq. 8. The ther-modynamic average is approximated as a time averageover 6 . × fm/c of simulation time. We presentresults for the angle averaged S ( q ) after averaging overapproximately 50 different directions of ~q . These re-sults are somewhat time consuming because we calculate S ( q ) for approximately 1400 different values of | ~q | for runrpcrust05.We calculate the inelastic contribution S ′ ( q ) by ap-plying a simple numerical filter that removes very rapidchanges in S ( q ) with q . Our filter, applied to a tableof q i and S ( q i ) values, works as follows: if S ( q i ) differsby more than some threshold ≈ . S ( q i − ) than g (r) Se (Z=34)O (Z=8)
FIG. 4: (Color on line) Radial distribution function g ( r ) ver-sus r over the mean ion sphere radius a , for run rpcrust-05 ata temperature T = 0 . OCPfit value is from the S ′ ( q ) fit in ref.[8] while Λ OCP is our calculation for simulation OCP.Γ Λ
OCPfit Λ OCP
250 0.362 0.348 q i and S ( q i ) are removed from the table. This removesnumerical noise and may remove delta function like con-tributions from the Bragg peaks. In addition, we maysimply miss some Bragg peaks because we only calculate S ( q ) for a finite number of q points. Our motivation forthis simple procedure is to calculate S ′ ( q ) and Λ basedon S ( q ) calculations that are not likely contaminated byBragg contributions.We first test this procedure with the one componentplasma simulation OCP of Table II, see Fig. 5. Our re-sults for S ′ ( q ) show more structure than the simple fitpresented in ref. [8]. Note that this may reflect a limita-tion of the fit. In addition, there is some high frequencynoise in our simulation. However, there is good agree-ment, to 4 %, between our OCP simulation and the fitfor the integral of S ′ ( q ) over q that is needed to calcu-late the Coulomb logarithm Λ, see Eq. 7 and Table III.Therefore, our procedure for S ′ ( q ) reproduces the knownCoulomb logarithm and thermal conductivity of a onecomponent plasma.Figures 6, 7, and 8 show S ( q ) for run rpcrust-05 at tem-peratures of T = 0 .
1, 0.2, and 0.3 MeV respectively. Weexpect similar results for run rpcrust-01b. These figuresalso show the simple fit to S ′ ( q ) for an OCP presentedin ref. [8]. This fit is significantly below S ′ ( q ) for run S ’( q ) OCPOCP fit
FIG. 5: (Color on line) Inelastic contributions to the staticstructure factor S ′ ( q ) versus momentum transfer q times themean ion sphere radius a , for run OCP, solid black line andthe simple fit presented in ref. [8] , dashed red line.TABLE IV: Coulomb Logarithm, Eq. 7, for rp process ashcomposition. The Λ values are from run rpcrust05 at theindicated temperatures T , while Λ OCPfit is from ref. [8] fora pure OCP and Λ
OCP+imp also includes impurity scatteringfrom ref. [11] with Q = 38 . T (MeV) Γ Λ OCPfit Λ OCP+imp
Λ0.1 850 0.104 0.146 0.1730.2 425 0.232 0.276 0.3660.3 283 0.334 0.377 0.530 rpcrust-05. Finally, these figures show the contributionof impurity scattering from [11] added to the OCP fit re-sults. Impurity scattering depends on Q , see Eq. 1 and Q = 38 . S ( q ) in Figs. 5-8 show statistical noise.However some of the effects of this noise average to zerowhen one integrates over S ′ ( q ) to calculate Λ. We es-timate the statistical error in our calculation of Λ at S ( q ) FIG. 6: (Color on line) Static structure factor S ( q ) for runrpcrust-05 versus momentum transfer q times the mean ionsphere radius a , dotted black line, at a temperature T = 0 . S ′ ( q ). This is calculated with a simple numericalfilter applied to S ( q ). Finally the dashed green line is the fitto OCP results for S ′ ( q ) from ref. [8] and the dashed dottedred line adds impurity scattering from ref. [11] to these OCPresults.TABLE V: Thermal conductivity for a temperature of T =0 .
043 MeV (5 × K). Results have been scaled to the in-dicated densities. The thermal conductivity κ is from runrpcrust-05 while κ OCP+imp is based on the OCP plus impu-rity values Λ
OCP+imp in Table IV and κ OCP is for the fit toan OCP without any impurity scattering.Γ ρ κ
OCP κ OCP+imp κ (g/cm ) (erg/K cm s) (erg/K cm s) (erg/K cm s)850 7 . × . × . × . ×
425 9 . × . × . × . ×
283 2 . × . × . × . × T = 0 . ± × fm/c to 4 × fm/c to a calculation usingconfigurations from 2 × to 3 × fm/c. We emphasizethat our procedure to calculate S ′ ( q ) from S ( q ) is modeldependent. Our numerical filter not only removes Braggpeaks but it may also remove some statistical noise. Notethat removing some noise seems to have minimal effectson the values of Λ that we calculate. We do not believeour results in Table IV are very sensitive to our procedureto determine S ′ ( q ). This is based on explicit calculationswith a few different procedures. S ( q ) FIG. 7: (Color on line) Static structure factor S ( q ) for runrpcrust-05 versus momentum transfer q times the mean ionsphere radius a , dotted black line, at a temperature T = 0 . S ′ ( q ). This is calculated with a simple numericalfilter applied to S ( q ). Finally the dashed green line is the fitto OCP results for S ′ ( q ) from ref. [8] and the dashed dottedred line adds impurity scattering from ref. [11] to these OCPresults. B. Thermal Conductivity
We now calculate the thermal conductivity κ using ourresults for the Coulomb logarithms. These results can bescaled to a range of densities n and temperatures T sothat the value of Γ, Eq. 3, remains the same. Table Vpresents κ at a temperature of T = 5 × K (a typicalvalue for a super bursting star). The thermal conduc-tivity is lower for run rpcrust05 than for an OCP. First,this is because run rpcrust05 has a large number of im-purities, corresponding to the large impurity parameter Q = 38 .
9. Second, we think κ may be further reduced be-cause the impurities in run rpcrust05 are not distributeduniformly. Instead they are concentrated in one region.Although our simulations show some of the effects ofimpurities on the thermal conductivity, we emphasizethat there may be important finite size effects becausewe find clustering. It is unrealistic to describe a largesystem by simply repeating our small simulation volumemany times. This would describe the impurities as beingconcentrated into many very small regions. Instead, webelieve the concentration of impurities indicates phaseseparation. We think that a large sample will separateinto two (or more) bulk phases. It is important to studyphase separation further with larger molecular dynamicssimulations and this may change our thermal conductiv-ity results. In general, one phase will be enriched in high Z ions while the other is enriched in low Z ions. Phase S ( q ) FIG. 8: (Color on line) Static structure factor S ( q ) versusmomentum transfer q times the mean ion sphere radius a forrun rpcrust-05, dotted black line, at a temperature T = 0 . S ′ ( q ). This is calculated with a simple numericalfilter applied to S ( q ). Finally the dashed green line is the fitto OCP results for S ′ ( q ) from ref. [8] and the dashed dottedred line adds impurity scattering from ref. [11] to these OCPresults. separation may act to reduce the impurity parameter Q and increase the thermal conductivity. For example, Q will be reduced in the high Z phase because low Z im-purities have gone into the other phase.In addition, nuclear reactions may reduce Q further.In general, we expect nuclear reactions to preferentiallyburn low Z impurities because of their low Coulomb bar-riers. See for example ref. [14]. This will reduce Q andincrease the thermal conductivity. One should study how Q evolves with depth because of reactions. Finally, it isimportant to analyze observations of crust cooling afterextended outbursts [3, 4] to see what observational con-straints can be placed on the thermal conductivity and Q . It may be that observations of rapid crust cooling canstrongly limit the size of Q [22].Phase separation may have another important effect.It will create layers in the crust of different compositionsand densities. These layers may not be spherically sym-metric. For example, phase separation depends on tem-perature. Therefore an anisotropic temperature distribu-tion will lead to an anisotropic density. It is importantto study how phase separation will change the structureof the star. IV. SUMMARY AND CONCLUSIONS
The crust of an accreting neutron star, likely, has acomplex composition with many impurities. Nuclei aresynthesized via the rapid proton capture process and thecomposition is modified by electron capture as material isburied to greater densities. We have performed MD sim-ulations, with a complex composition, to study the struc-ture of the crust. Our simulations form ordered crystalsrather than an amorphous solid.However, we find phase separation. Some low Z impu-rities are concentrated into a subregion of the full simu-lation volume. This phase separation, between two solidphases, is similar to the chemical separation found pre-viously between liquid and solid phases [12]. Previously,we assumed a steady state equilibrium where chemicalseparation greatly increases the concentration of low Z impurities in the liquid ocean. However, the compositionof the solid crust was assumed to be the same as thatof the accreting material. Our new results disprove thissteady state assumption.The crust can not be uniform, given our initial com-position. Phase separation will divide the crust into twoor more regions of different compositions. This may haveimportant implications for the structure of the star. Forexample, composition anisotropies could lead to gravi-tational wave radiation from a quadrupole deformation [23]. In future work we will study the size of possible com-positional asymmetries because of an anisotropic temper-ature distributionWe calculated the static structure factor S ( q ) for oursimulations and from S ( q ) the thermal conductivity κ .Since our simulations have a complex composition, weautomatically include the contributions of impurity scat-tering. We find that κ is somewhat reduced because ofimpurity scattering and because the impurities are notdistributed uniformly. We expect the same results forthe electrical conductivity σ , that is important for mag-netic field decay [26], and the shear viscosity η , that candamp neutron star oscillations [24],[25]. The reductionin κ may be observable in crust cooling times and theseobservations may set limits on impurities. Future workshould study how phase separation and or nuclear reac-tions impact impurity concentrations. V. ACKNOWLEDGMENTS
We thank Ed Brown and Andrew Cumming for helpfuldiscussions and acknowledge the hospitality of the Insti-tute for Nuclear Theory where this work was started.This work was supported in part by DOE grant DE-FG02-87ER40365 and by Shared University Researchgrants from IBM, Inc. to Indiana University. [1] R. Wijnands et al., astro-ph/0405089.[2] E. M. Cackett et al., MNRAS , 479 (2006).[3] R. E. Rutledge et al., ApJ. , 413 (2002).[4] P. S. Shternin, D. G. Yakovlev, P. Haensel, and A. Y.Potekhin, MNRAS , L43 (2007).[5] A. Cumming and L. Bildsten, ApJ (2001) L127.[6] T. E. Strohmayer and E. F. Brown, ApJ (2002) 1045.[7] A. Cumming, J. Macbeth, J. J. M. in ’t Zand and D.Page, ApJ. (2006) 429.[8] A. Y. Potekhin, D. A. Baiko, P. Haensel, and D. G.Yakovlev, Astron. Astrophys., (1999) 345.[9] H. Schatz et al., PRL (2001) 3471.[10] S. Gupta, E. F. Brown, H. Schatz, P. Moller, and K-L.Kratz, ApJ (2007) 1188.[11] N. Itoh and Y. Kohyama, ApJ. (1993) 268.[12] C. J. Horowitz, D. K. Berry, and E. F. Brown, PRE (2007) 066101.[13] E. F. Brown, private communication.[14] C. J. Horowitz, H. Dussan and D. K. Berry, PRC (2008) 045807.[15] S. E. Woosley, A. Hager, A. Cumming, R. D. Hoffman, J.Pruet, T. Rauscher, J. L. Fisker, H. Schatz, B. A. Brown,and M. Wiescher, ApJ Supp. (2004) 75.[16] H. E. Dewitt, W. L. Slattery, and J. Yang in “StronglyCoupled Plasmas”, eds. H. M. Van Horn and S. Ichimaru,Univ. of Rochester Press 1993, p425.[17] D. G. Yakovlev, L. R. Gasques, M. Wiescher, and A. V.Afanasjev, PRC (2006) 035803.[18] B. Jancovici, J. Stat. Phys. (1977) 357.[19] N. Itoh, S. Uchida, Y. Sakamoto, and Y. Kohyama, arXiv:0708.2967.[20] L. Verlet, Phys. Rev. (2000) 902.[24] A. I. Chugunov and D. G. Yakovlev, Astronomy Reports (2005) 724.[25] C. J. horowitz and D. K. Berry, PRC (2008) 035806.[26] P. Goldreich and A. Reisenegger, ApJ.395