Relaxation by thermal conduction of a magnetically confined mountain on an accreting neutron star
Arthur G Suvorov Andrew Melatos




firstpage–lastpage ?

Possibility of negative refraction for visible light in disordered all-dielectric resonant high index metasurfacesalink
    Andrey V. Panov
    Accepted ?. Received ?; in original form ?
=======================================================================================================================





A magnetically confined mountain on the surface of an accreting neutron star simultaneously reduces the global magnetic dipole moment through magnetic burial and generates a mass quadrupole moment, which emits gravitational radiation. Previous mountain models have been calculated for idealized isothermal and adiabatic equations of state. Here these models are generalised to include non-zero, finite thermal conduction. Grad-Shafranov equilibria for three representative, polytropic equations of state are evolved over many conduction time-scales with the magnetohydrodynamic solver PLUTO. It is shown that conduction facilitates the flow of matter towards the pole. Consequently the buried magnetic field is partially resurrected starting from an initially polytropic Grad-Shafranov equilibrium. The poleward mass current makes the star more prolate, marginally increasing its detectability as a gravitational wave source, though to an extent which is likely to be subordinate to other mountain physics. Thermal currents also generate filamentary hot spots (≳ 10^8 K) in the mountain, especially near the pole where the heat flux is largest, with implications for type I X-ray bursts. 




stars: neutron, accretion, magnetic fields, gravitational waves




§ INTRODUCTION
 

Observations of binary neutron stars with white dwarf or supergiant companions and a history of accretion suggest that the neutron star magnetic dipole moment μ decreases over time, as the accreted mass M_a increases <cit.>. Several theoretical mechanisms exist to explain the trend, such as accelerated Ohmic decay <cit.>, interactions between superfluid vortices and superconductor flux tubes within the stellar interior <cit.>, or the process of magnetic burial <cit.>. In magnetic burial, the focus of this paper, matter is guided onto the polar cap by the magnetic field to form a mountain-like density profile supported by the compressed equatorial magnetic field (`magnetic mountain') <cit.>. The resulting mass quadrupole moment emits gravitational radiation <cit.>.










The short-term stability and long-term relaxation of a magnetic mountain have been studied by several authors. In the short term, on the Alfvén and tearing-mode time-scales, axisymmetric mountain equilibria are susceptible to the undular submode of the Parker instability <cit.> and to pressure-driven toroidal-mode instabilities <cit.>, once M_a exceeds a critical threshold. The system is not necessarily disrupted; the instability saturates, and the mountain adjusts to a new equilibrium, stabilized by magnetic line-tying at the stellar surface and the compressed magnetic `wall' at the equator <cit.>. In the long term, the mountain relaxes due to Ohmic dissipation <cit.>, soft-crust sinking <cit.>, or a combination of the latter two processes <cit.>. Its structure is modified also by factors like the Hall effect <cit.> and the equation of state (EOS) <cit.>. 

Mountains on recycled pulsars may be responsible for the discrepancy between magnetic field strengths inferred from spin-down and cyclotron line measurements <cit.>. The local magnetic field can be ∼ 10^4 times stronger than the global value inferred from μ <cit.>. Once M_a increases beyond a certain level, phase-dependent cyclotron resonance scattering features are predicted to emerge in the X-ray spectrum <cit.>. Additionally, X-ray observations of neutron star binaries reveal type I X-ray bursts with recurrence times ranging between a few minutes and ∼ 10^3 hours <cit.>. Recurrence times ≲ 10 min (e.g. in 4U 1608–522) are too short for many theoretical ignition models and may indicate the existence of multiple, isolated patches of fuel on the stellar surface <cit.>, which are fenced-off magnetically if the polar magnetic field geometry is complicated <cit.>. A simultaneous detection of gravitational waves, X-ray bursts with short recurrence times, and cyclotron features in an X-ray binary some time in the future would strongly indicate the presence of a magnetic mountain <cit.>. 


In this paper we include thermal conduction in magnetic mountain models self-consistently for the first time. Thermal conduction is potentially important, because the mountain forms at an elevated temperature, caused by accretion-driven heating, and cools through its sides (if accretion is confined to a narrow column) or throughout its volume (once accretion switches off). Thermal fluxes directed out of localized polar hot spots control the instantaneous hydromagnetic structure of the mountain by regulating the EOS <cit.> and the long-term, quasistatic relaxation of the mountain by regulating temperature-sensitive dissipative mechanisms like Ohmic decay <cit.>. Modeling thermal conduction in magnetic mountains self-consistently is therefore important for understanding the relationship between hot spots, magnetic fields, X-ray burst activity, and gravitational radiation, providing the basis for multi-messenger tests of the polar magnetic burial scenario. 


The purpose of this paper is to elucidate, with the aid of numerical simulations, the dominant thermal processes that modify the short-term structure and long-term evolution of a magnetic mountain, when thermal conduction is “switched on” in the model. Predictions are made, in broad qualitative terms, regarding how potentially observable properties (e.g. μ) are affected by thermal conduction. We emphasize, however, that the simulations are not yet at the point where they yield highly realistic mountain models, which are ready to be compared in detail with observational data. Such comparisons would require a more sophisticated description of the stratified structure, composition, and EOS of the crust, better observational knowledge of the high-order magnetic multipoles near the surface, and expanded computational resources to handle the disparate thermal and hydromagnetic time-scales in the problem. Our investigation proceeds in two stages. In Section 2, we use the Grad-Shafranov solver developed by <cit.> and extended by <cit.> to calculate the steady-state structure given an EOS and an initially dipolar magnetic field. We then numerically evolve the Grad-Shafranov equilibrium using the magnetohydrodynamics (MHD) code  <cit.> with and without thermal conduction in Sections 3 (set-up details and local mountain properties) and 4 (global observables) and compare the effects on potentially observable properties such as μ. Long-term thermal relaxation is explored in Section 5. Finally, the astrophysical implications of the results, including for gravitational wave emission, are discussed briefly in Section 6.










§ POLAR MAGNETIC BURIAL




 §.§ Qualitative behaviour


During accretion, the neutron star's polar magnetic field buckles underneath the infalling matter, and the field lines spread equatorially due to flux freezing. The lateral pressure gradient at the base of the accreted mountain is balanced by the Lorentz force in the compressed, equatorial magnetic belt. This process is illustrated schematically in Figure <ref> [see also Figure 6 of <cit.>]. The compressed magnetic field is more intense than the pre-accretion field locally, due to magnetic flux conservation, but the global moment μ reduces, because the magnetic distortion induces screening currents, which reduce the radial magnetic field near the pole <cit.>. 

It is observed that μ decreases with M_a in binary systems <cit.>. <cit.> proposed the widely used, empirical law 

    μ = μ_i( 1 + M_a / M_c)^-1 .

In (<ref>), we define M_c to be the critical accreted mass, for which the global dipole moment is halved. The dipole moment before accretion begins is given by μ_i = B_⋆ R_⋆^3. We take R_⋆ = 10^6 cm for the stellar radius and B_⋆ = 10^12.5 G for the natal magnetic field strength at the polar surface, in line with population synthesis models <cit.>. Self-consistent, MHD simulations reproduce the empirical scaling (<ref>) for small accreted masses M_a≪ M_c in isothermal and adiabatic mountains with M_c∼ 10^-6 M_⊙ and M_c∼ 10^-8 M_⊙ respectively <cit.>. For 10^-1≲ M_a/ M_c≲ 10, the simple estimate in (<ref>) breaks down and the burial effect is better represented by a power-law μ / μ_i = (M_a / M_c)^-a, where 1 ≤ a ≤ 2.47  depends on the EOS [see section 4.1 of <cit.> and Fig. 8(c) of <cit.>]. Numerical difficulties prevent simulations from probing the regime M_a / M_c≳ 10, where a significant deviation from (<ref>) is expected <cit.>, though Ohmic diffusion sets a burial limit of μ / μ_i≳ 10^-8 <cit.>.






The critical mass M_c depends strongly on the EOS <cit.>. For a softer EOS, the mountain has a relatively small thickness (∼ 10^3 cm), because the material is easier to compress. Strong local magnetic fields (≲ 10^15 G) exist near the stellar surface, as polar field lines buckle, and the polar magnetic flux is squeezed into a small volume. Consequently, the screening currents flow closer to the stellar surface for a softer EOS than for a harder EOS, and μ reduces less for a given M_a. For a softer EOS, the critical mass is found to lie in the range 10^-5≲ M_c / M_⋆≲ 10^-2 <cit.>. For a harder EOS, the mountain is thicker (∼ 10^4 cm), and μ reduces further. For a polytropic EOS with index Γ≳ 4/3, the critical mass is found to lie in the range 10^-9≲ M_c / M_⋆≲ 10^-6 <cit.>. <cit.> found M_c∝ B_⋆^2. 











 






 
 
The mountain mass quadrupole moment can be expressed in terms of the mass ellipticity  ϵ, which is given approximately by <cit.>

    ϵ≈( 5 M_a / 4 M_⋆) ( 1 + 9 M_a / 8 M_c)^-1 ,





where M_⋆ = 1.4 M_⊙ is the stellar mass. Therefore, for a given accreted mass, there is a one-to-one relationship between μ and ϵ through (<ref>) and (<ref>) which depends on the value of M_c. The mass quadrupole moment emits gravitational radiation, as the star spins <cit.>. The implications are discussed in Section 4.2.








 §.§ Hydromagnetic equilibrium


The hydromagnetic structure of the mountain has been calculated for various EOS previously <cit.>. In this paper, we start by solving for the structure of a steady-state (∂ / ∂ t = 0) and immobile (= 0, where v is the fluid velocity) mountain. We then input the result into  as the starting point for time-dependent simulations which include thermal conduction. In the first stage we do not model growth of the mountain but rather solve for a self-consistent equilibrium (i.e. hydromagnetic force balance) given a certain amount of accreted mass and an EOS. Time-dependent simulations in the literature confirm that the equilibrium agrees closely with mountains built from scratch by injecting mass from below <cit.>.

We assume that the magnetic field  may be described through an axisymmetric Chandrasekhar decomposition without a toroidal component[ Three-dimensional mountain equilibria for M_a≳ M_c are susceptible to Parker-like <cit.> and ballooning <cit.> instabilities with EOS-dependent growth rates <cit.>. However, the instabilities do not disrupt the mountain; they reduce the ellipticity by ≲ 30 percent when they saturate <cit.>. Furthermore, previous three-dimensional, time-dependent simulations reveal that the magnetic field relaxes to an almost axisymmetric configuration within a few Alfvén times <cit.>.] for simplicity, i.e. B takes the form

    = ψ/r sinθ×ê_φ

in spherical coordinates (r,θ,φ), where ψ(r,θ) is a scalar flux function <cit.>. The mountain equilibrium is determined by the force balance (Grad-Shafranov) equation <cit.>

    p + ρϕ + Δ^2 ψψ = 0,

where Δ^2 denotes the Grad-Shafranov operator,

    Δ^2 = 1/4 π r^2 sin^2θ[ ∂^2/∂ r^2 + sinθ/r^2∂/∂θ( 1/sinθ∂/∂θ) ],
 
p is the fluid pressure, ρ is the mass density, and ϕ is the gravitational potential. For now, we assume that the accreted matter obeys a barotropic EOS, p = p(ρ), once equilibrium is reached. Although neutron stars are expected to be non-barotropic <cit.>, the barotropic and non-barotropic solutions to the Grad-Shafranov problem are broadly similar, even in magnetars <cit.>.



For p = p(ρ), equation (<ref>) implies that ρ and ψ satisfy

    0 = d p/d ρ( ∂ρ/∂ r∂ψ/∂θ - ∂ρ/∂θ∂ψ/∂ r) + ρ( ∂ϕ/∂ r∂ψ/∂θ - ∂ϕ/∂θ∂ψ/∂ r).

Equation (<ref>) can be solved exactly using the Lagrange-Charpit method to yield <cit.>

    ∫d p/d ρd ρ/ρ = F(ψ) - ( ϕ - ϕ_0),

where ϕ_0 denotes a reference gravitational potential at the neutron star surface, and F is an arbitrary function of the scalar flux. Given p = p(ρ), the integral in (<ref>) can be evaluated, in principle, to express ρ in terms of ψ or vice-versa. 

In order to obtain a one-to-one correspondence between the pre- and post-accretion states that respects flux-freezing, we demand that the steady-state, mass-flux ratio d M / d ψ, defined as the mass enclosed between the infinitesimally separated flux surfaces ψ and ψ + d ψ, equals that of the initial state plus any accreted matter <cit.>.  This restriction on M(ψ) leads to the constraint <cit.>

    d M/d ψ = 2 π∫_C ds ρ[ r(s), θ(s) ] r sinθ | ψ |^-1 ,

where C is the curve ψ[r(s),θ(s)] = ψ parametrized by the arc length s. Equation (<ref>) can be solved by inverting (<ref>) for ρ in terms of ψ given p=p(ρ), leading to a unique expression for the function F(ψ) given M(ψ). The explicit forms of F for adiabatic and isothermal EOS can be found in <cit.> and <cit.> respectively. 



We solve (<ref>) simultaneously with (<ref>) numerically using the relaxation algorithm described in <cit.> and later extended by <cit.>. Specifically, the solver employs iterative under-relaxation combined with a finite-difference Poisson solver to solve (<ref>) for ψ, obtain ρ from (<ref>), and feed the result back into (<ref>) iteratively, until convergence is achieved. Additional information regarding units, convergence, and stability can be found in the aforementioned papers and is not repeated here; see also <cit.> and <cit.>. In accord with previous work, we prescribe the mass-flux distribution in one hemisphere to be <cit.>

    M(ψ) = M_a[ 1 - exp(-ψ/ψ_a)]/2 [ 1 - exp(-b)],

where M_a is the accreted mass, ψ_a labels the field line at the polar-cap boundary (that closes just inside the inner edge of the accretion disc), and we define b = ψ_⋆ / ψ_a, where ψ_⋆ labels the total hemispheric flux. Throughout this paper we set b=3 to ensure numerical stability.



For simplicity we assume a constant gravitational acceleration, with ϕ(r) = G M_⋆ r / R_⋆^2, where R_⋆ is the stellar radius, and make the Cowling approximation (i.e. we ignore self-gravity). These assumptions are justified, because the mountain never rises more than ∼ 10^4 cm above the surface at r = R_⋆, and we consider systems with M_a / M_⊙≲ 10^-1 [see section 2.1 of <cit.>]. Additionally, <cit.> and <cit.> found that the Cowling approximation alters the mass ellipticity by at most a factor of ∼ 3 even for the strongest magnetar fields (≲ 10^16 G).





Equation (<ref>) is solved subject to physically motivated boundary conditions, which carry through to the evolution experiments in  in Sec. 2.4 (see also Appendix A). Following previous work, we set ψ(R_in,θ) = ψ_⋆sin^2θ (surface dipole), ∂ψ / ∂ r (R_m,θ) = 0 (Neumann outflow[Ideally, this condition would be replaced by a dipolar field condition at the outer edge, i.e. ψ(R_m,θ) = ψ_msin^2θ for some value of ψ_m, so as not to introduce artificial magnetic multipoles (including a monopole) far from the stellar surface. However, in order to assign a value to ψ_m, as necessary for numerical computation, we need to know by how much magnetic burial reduces the dipole moment, with μ / μ_i∝ψ_m / ψ_i. In principle, it is possible to adjust ψ_m iteratively in order to obtain a self-consistent simulation, but this is technically challenging (see also Footnote 4). A thorough discussion of the issue can be found in Sec. 4.3 of <cit.> and in Sec. 4.1 of <cit.>, as well as references therein. In particular, <cit.> showed that the density distribution (Fig. 14 of the latter reference, left panel) is virtually indistinguishable between Neumann and dipole condition equilibria, while the magnetic field lines (right panel) tend to agree except in the outermost regions, where the plasma density is low. The Neumann outflow condition artificially increases μ, relative to corresponding  equilibria, by ≲ 10 % [see Fig. 3(f) of <cit.>].]), ψ(r,0) = 0 (straight polar field line), and ∂ψ / ∂θ (r , π/2 ) = 0 (equatorial symmetry), where R_in≤ r ≤ R_m and 0 ≤θ≤π/2 demarcate the computational volume <cit.>. The outer radius R_m is chosen large enough to encompass all the screening currents and the outer edge of the accreted matter; we set R_m = 1.4 R_⋆ throughout this paper. The boundary conditions on ψ are reformulated as conditions on  through (<ref>) and conditions on ρ (and hence p) through (<ref>). In the time-dependent  simulations (see Sec. 2.4) we also stipulate no slip at R_in, outflow at R_m, and reflecting boundary conditions on v at the equator.

Although R_in is treated as a hard surface for simplicity, it is not so in reality; a mountain several tens of meters high, whose base reaches neutron drip densities, sinks into the lower-density substrate <cit.>. A full treatment of sinking requires time-dependent simulations. <cit.> showed that the results are approximated reasonably by hard-surface solutions (the mountain ellipticity decreases by a factor ≲ 2 for soft crust solutions relative to hard-surface solutions), if R_in corresponds to the layer above which the mass equals M_a; i.e. one has R_in < R_⋆ and M_a = ∫^r=R_⋆_r=R_in d^3x ρ. As R_in is fixed, the stellar radius R_⋆ varies slightly (≤ 1 per cent for M_a≤ 10^-4 M_⊙) between models with different M_a but the same EOS. Figure <ref> visually demonstrates the relationship between R_in, R_⋆, and M_a. 









 §.§ Equation of state


Crustal matter experiencing compression due to accretion undergoes a variety of non-equilibrium nuclear processes, such as electron captures and beta decay, neutron emission and absorption, and pyconuclear fusion, all of which play a role in determining the EOS of the accreted crust <cit.>. The original outer crust, consisting of cold, catalysed matter, is replaced by a new, non-catalysed crust after ≲ 10^5 yr <cit.>. The EOS of an accreted, non-catalysed crust, relevant for our calculations, has been calculated numerically by <cit.> by modeling the non-equilibrium processes listed above, using the compressible liquid drop model of <cit.> to estimate the various thermodynamic rates which feed into the Gibbs equation [see also Sec. 2.4 of this paper and <cit.>].




In this paper we consider three idealized yet physically motivated polytropic EOS: a single-index model which best approximates (in a least-squares sense; see below) a realistic accreted crust (model A), and two of the classical ideal gas models considered by <cit.>, corresponding to a gas of non-relativistic degenerate electrons (model B), and a gas of non-relativistic degenerate neutrons (model C). Their parameters are summarised in Table <ref>; see also <cit.>. Models A, B, and C apply only at t=0; they are used to construct initial conditions for . The Grad-Shafranov solver developed by <cit.> quasi-statically determines the `end-state' of an adiabatic accretion process in the absence of thermal conduction. As noted above, we initialise  with an M_a-dependent Grad-Shafranov equilibrium to avoid numerical difficulties, cf. <cit.>. In reality, however, the true end-state of accretion depends on Ṁ_a as well as M_a. Accreted plasma on the stellar surface is expected to be approximately isothermal (Γ≈ 1) at all depths for low-accretion rates Ṁ_a≲ 10^-10 M_⊙ yr^-1 <cit.>. In contrast, a crust formed on a star accreting near the Eddington limit Ṁ_a≳ 10^-8 M_⊙ yr^-1 has a more complicated polytropic EOS with a depth-dependent adiabatic index (1 ≲Γ≲ 5/3) <cit.>. Hence the time-dependent accretion process, which is not modeled here except implicitly through (<ref>), directly affects the softness or hardness of the EOS. In particular, a self-consistent model of accretion along the lines of the sinking problem treated by <cit.> would lead to end-state values of k_Γ and Γ which depend on both M_a and Ṁ_a <cit.>. The numerical experiments we conduct in , which evolve the EOS via thermal conduction, partially account for the effects of a near-Eddington accretion rate on the crustal EOS a posteriori (see Sec. 2.4). This procedure has been validated in the absence of thermal conduction by <cit.>. 


In Figure <ref> we graph pressure-density relationships for models A, B, and C (broken curves) together with the numerical results of <cit.> (solid curve). For ρ≳ 1.5 × 10^13 g cm^-3, the maximum density computed by <cit.>, we graph the inner-crust model of <cit.>, also computed using the compressible liquid drop model <cit.>. Denoting the neutron drip density by ρ_nd [ρ_nd∼ 5 × 10^11 g cm^-3 in an accreted crust <cit.>], we see that the numerical results are approximated adequately by models B and C used in previous work <cit.> in the regimes ρ≪ρ_nd and ρ≳ρ_nd respectively. On the other hand, model A is constructed to uniformly approximate the realistic EOS for all ρ≲ 10^14 g cm^-3. The parameters k_Γ and Γ for model A are computed by fitting p(ρ) = k_Γρ^Γ with the Levenberg-Marquardt (damped least-squares) algorithm <cit.> to the data collated in Table 1 of <cit.>. Denoting the <cit.> numerical pressure by p^HZ and the model A pressure by p^A, the fit yields relative errors of 0.94 ≤ p^HZ/p^A≤ 2.51 for 10^8≤ρ/  g cm^-3≤ 10^14. Throughout most of the mountain volume by mass, i.e. for ρ≥ 10^12 g cm^-3, the errors drop to ≤ 6 percent, with 0.95 ≤ p^HZ/p^A≤ 1.06.  For ρ≤ 10^9 g cm^-3 we have 0.94 ≤ p^HZ/p^A≤ 1.07.












A piecewise polytropic fit to the solid curve in Fig. <ref> (e.g. a spline fit to log p vs logρ) is a better approximation than the uniform, Γ = 1.18 fit in model A. As a practical matter, however, it is difficult to generalise the Grad-Shafranov calculation in Sec. 2.2, especially the explicit formula for F(ψ) [equation (8) in <cit.>], to apply to multiple layers with proper matching at the inter-layer boundaries. As the Grad-Shafranov calculation is an essential input to the  simulations, the uniform Γ = 1.18 approximation is the best we can do for now. For this reason, among others, the final results should be viewed as qualitatively representative models of the thermal conduction physics rather than quantitatively accurate mountain models to be compared in detail to observational data.

Different EOS predict different maximum (base) densities ρ_max and heights H for any given M_a. In Table <ref> we list the characteristic ρ_max and H for runs performed in this paper (see Sec. 3.2) together with a rough estimate for the expected depth within a neutron star to which ρ_max corresponds. Note that the sinking depths listed in Table <ref> apply for stellar mass M_⋆ = 1.4 M_⊙, and a Skyrme-type EOS at zero temperature, used to describe both the crust and the liquid core, based on the effective nuclear interaction SLy <cit.>. Different EOS and stellar masses lead to different sinking depths. For models A and C, the theoretical depth exceeds the simulated height of the mountain. Depending on the crustal elasticity <cit.>, this indicates that mountain matter should sink beneath the surface and influence the hydromagnetic structure of the star <cit.>. However, using breaking strain arguments, it has been shown that realistic crustal ellipticites of neutron stars cannot exceed ∼ 10^-6 <cit.>, which is less than those associated with accreted mountains (see Sec. 4.2). As such, any gravitational radiation due to crustal quadrupole moment generation via back-reaction effects from a sinking mountain is likely to be dwarfed by the radiation due to the mountain itself <cit.>, though there may be interesting consequences for other phenomena, e.g. crust-core coupling <cit.>. The lateral (θ) structure of the mountain is not affected greatly by sinking, as shown by <cit.>. Hence the main effect of sinking on ϵ is to reduce it by a factor ≈(R_pre / R_post)^5, where R_pre and R_post are the characteristic radii of the base of the mountain before and after sinking, respectively, and “before sinking" here means “in the context of a hard-surface Grad-Shafranov calculation”. In any case, because we do not model sinking, the values of the ellipticities (and heights) presented in this paper should be taken as upper limits. Modeling a realistic neutron star together with a sinking mountain in a way that simultaneously tracks the Alfvén and sinking time-scales is a difficult problem that will be considered in future work.







 §.§ MHD evolution

The steady-state solution to the Grad-Shafranov problem in section 2.2 serves as initial data for evolving the mountain dynamically. In the absence of viscosity and under the assumptions of infinite electric conductivity[The Ohmic diffusion and thermal conduction time-scales (see Sec. 2.5) are in the ratio τ_d / ≈ 2 × 10^-14( σ /  s^-1) || / |^2  | L^-2( ρ / 10^10 g cm^-3)^-1, for characteristic length-scale L and electrical conductivity σ. In the crust-magnetosphere interface, one has σ≲ 10^16 s^-1 <cit.>. In the inner crust one has σ≳ 10^24 s^-1 <cit.>. Hence we find τ_d / ≫ 1 throughout the computational volume for the range of accreted masses considered in this paper, even in regions with strong magnetic gradients, because the density is low there (ρ≲ 10^10 g cm^-3). We can therefore safely ignore the effects of Ohmic diffusion over the time-scales simulated within this paper; see also <cit.>.] (ideal MHD) and the Cowling approximation, the evolution is governed by the continuity, Euler, and Faraday equations, which read <cit.>

    0 = ∂ρ/∂ t + ·( ρ),


    0 =   ρ∂/∂ t + ρ( ·)  +  p + ρϕ
        -   ( 4 π)^-1[ (×) ×] ,

and

    0 = ∂/∂ t - ×( ×),




respectively, given ϕ. The MHD equations (<ref>)–(<ref>) are closed by the energy equation,

    ·F =   ·[ ( ε + ρ^2/2 + p + ^2/4 π + ρϕ)   - ( ·)/4 π] 
       +∂/∂ t( ε + ρ^2/2 + ^2/8 π + ρϕ),

where T is the temperature, F is the heat flux, and ε is the internal energy <cit.>. The kinetic properties of the fluid determine the internal energy in terms of the thermodynamic variables p, ρ, and temperature T, i.e. ε = ε(p,ρ,T) <cit.>. The Gibbs fundamental equation,

    0    = T d S - p d V - d ε
       = T d S - p d V - ( ∂ε/∂ p dp + ∂ε/∂ρ dρ + ∂ε/∂ T dT ),

where V is the system volume and S is the entropy, provides an additional constraint for the state variables. We then have seven scalar equations, namely (<ref>)–(<ref>), for seven variables in the axisymmetric problem: p, ρ, T, B_r, B_θ, v_r, and v_θ. In practice, (<ref>) determines T given p and ρ, while (<ref>) determines the the relationship between p and ρ, i.e. the barotropic EOS p = p(ρ) initially. Under the assumption of an ideal gas (consistent with a polytropic EOS), equation (<ref>) leads to the well-known relation <cit.>

    T = η m_u p/k_Bρ,

where η is the mean molecular weight, m_u is the atomic mass unit and k_B is the Boltzmann constant. The temperature field T(p,ρ) in (<ref>) is calculated from the Grad-Shafranov output and forms an additional input into  for thermal conduction simulations. We make the assumption of symmetric nuclear matter to determine η for simplicity as in previous work [see section 2.3 of <cit.>]. Although the ideal gas law (<ref>) is modified in degenerate matter, it provides a good approximation for partially-degenerate, accreted material on a neutron star crust <cit.> and is straightforward to handle within . A sensitivity analysis associated with expression (<ref>) is presented in Appendix B, where it is shown that using (<ref>) to determine T as opposed to a degenerate EOS calculated from first principles overestimates the temperature by ≲ 15% throughout the bulk of the mountain (see also Secs. 2.3 and 3), where thermal transport matters most; a small effect compared to other uncertainties in the problem.









 §.§ Thermal conduction





In the presence of thermal conduction, the flux on the left-hand side of (<ref>) takes the form <cit.>

    F = κ_|| ( · T ) + κ_⊥ [  T -  ( · T ) ] ,

where the thermal conductivities κ_|| and κ_⊥, both measured in units of erg s^-1 K^-1 cm^-1, describe heat transport parallel and perpendicular to the magnetic field respectively. The conductivities of a magnetised, fully ionized plasma are dominated by electron transport. They are given in the diffusion approximation by the Balescu-Braginskii formulas [see <cit.> and Table 3.2 of <cit.>],


    κ_||≈  1.77 × 10^16( T_e/10^9  K)^5/2erg s^-1 K^-1 cm^-1 ,

 and

    κ_⊥≈    κ_||[ 1 + 1.19 ( ρ/10^10 g cm^-3)^-2( ||/5 × 10^13 G)^2
       ×( T_e/10^9  K)^3]^-1erg s^-1 K^-1 cm^-1 ,

assuming the Coulomb logarithm satisfies lnΛ = 30 ≈constant <cit.>. The estimates (<ref>) and (<ref>) are independent of the ideal gas assumption (<ref>). In the limit of a vanishing magnetic field, one has κ_⊥ = κ_||, and (<ref>) reduces to F = κ_|| T.
















 §.§ Time-scales

A mountain with the structure in Fig. <ref> contains steep density and magnetic field gradients, so there is no unique definition for characteristic time-scales, like the Alfvén time  and thermal conduction time . In order to analyse our numerical results in Sec. 3 onwards, we adopt the definition <cit.> 

    = L √(4 πρ̅)/ || 
    
       ≈ 7 × 10^-5( L/10^4  cm) ( ρ̅/10^10 g cm^-3)^1/2
             ×( ||/5 × 10^13 G)^-1 s,

where L is taken to be the density scale-height where ρ drops to 10^-3 times its maximum value ρ_max, ρ̅ is the volume-averaged density

    ρ̅ = 1/∫ d^3x ρ,

where  is the mountain volume (10^-7ρ_max≤ρ≤ρ_max) and || is the volume-averaged magnetic field strength,

    || = 1/∫ d^3x  ||.

All the quantities (<ref>)–(<ref>) are computed at t=0 to define  for any given run. Similarly, for thermal conduction, from the heat equation we have

    =   5/3( L / 1 cm)^2 ( ρ̅ /  g cm^-3) 
       ×[ κ̅_||  / erg s^-1 K^-1 cm^-1]^-1 s
    
       ≈ 94.2 ( L/10^4  cm)^2( ρ̅/10^10 g cm^-3) ( T̅_e/10^9  K)^-5/2 s,

where κ̅_⊥, κ̅_||, and T̅_e are volume-averaged quantities calculated in the same manner as ρ̅ in (<ref>) and || in (<ref>). 

The thermal conduction time is ∼ 10^6 times longer than the Alfvén time for a typical, realistic mountain. Computational expense restricts us to t ≲ 2 throughout most of this paper (see Sec. 3). The ratio of the time-scales varies from run to run and for different EOS. We perform some `long-term' evolutions (up to ∼ 80) in Section 5 to explore the effects of thermal relaxation. In Sections 3 and 4 we show that (<ref>) and (<ref>) agree with the time-scales of characteristic behaviours observed empirically in the simulations.





§ THERMAL EVOLUTION


PLUTO <cit.> is a general-purpose MHD solver designed to handle steep gradients associated with strong shock phenomena in astrophysical applications. It solves  (<ref>)–(<ref>)  given the Grad-Shafranov solution and the boundary conditions described in Sec. 2.2 as inputs. The details of the computation are presented in Appendix A along with convergence tests.

In this section we present results from   simulations of magnetic mountain evolution on time-scales comparable to .  We load a Grad-Shafranov equilibrium calculated in Sec. 2.2 for some equation of state (e.g. A,B, or C in Table <ref>) into  and evolve it in two ways, with thermal conduction switched on or off, to explore the mass density and magnetic field profiles for a variety of runs (Sec. 3.2), the evolution of the thermal flux (Sec. 3.3), and the evolution of global observables like μ (Sec. 4.1) and ϵ (Sec. 4.2). 



Several numerical and physical issues affect the  output. (i) The Grad-Shafranov code computes ψ, while  accepts the components of B. The calculation of B from ψ involves differentiation, which introduces some numerical error. We use 's inbuilt bi-linear interpolation algorithm to map the Grad-Shafranov output to  input (see Appendix A). (ii)  maintains ideal-MHD flux freezing through a Godunov scheme [e.g. <cit.>], but it does not act directly to satisfy the integral constraint (<ref>) on M(ψ). As M(ψ) is not an input into  , and equation (<ref>) is a non-linear equation for ψ, it is possible that multiple, valid solutions for ψ exist at any given t. One can imagine  picking a solution branch unpredictably based on numerical fluctuations, if two valid solutions for ψ are numerically close. (iii) The Grad-Shafranov equilibrium may not represent the stable endpoint of a well-posed initial value problem because the Grad-Shafranov equation has multiple unstable solutions <cit.>. This is related to the loss-of-equilibrium phenomenon investigated by <cit.>. 







 §.§ Representative example


We start by considering a representative simulation, which demonstrates the main features of thermal evolution: EOS model A with M_a = 1.8 × 10^-5 M_⊙≈ 0.58 M_c. We set up a polytropic initial state (Γ = 1.18 at t=0), allow thermal conduction to take place, and evolve the mountain. The EOS parameters, described in Table <ref>, are entered into the Grad-Shafranov solver, which produces the initial input for . The thermal profile is entered according to (<ref>). Two separate   instances are evolved, with and without the thermal flux F appearing in the right-hand side of equation (<ref>). The time-scales (<ref>) and (<ref>) read ≈ 2.0 × 10^-5 s and ≈ 5.6 × 10^6. 



In Figure <ref>, which demonstrates several features typical of an initially polytropic mountain, we graph contours of ρ and magnetic field lines at t=0. The mountain reaches a maximum altitude of ≈ 1.3 × 10^3 cm near θ = π/4 (where || rises to a maximum). Most of the mass (≈ 92 %) is concentrated within the octant 0 ≤θ≤π/4. The densest point, with ρ = 2.2 × 10^13 g cm^-3, lies at the pole. The magnetic field lines are shifted equatorially; the maximum contour lies at θ = π/4, in contrast to the initial dipole field (maximum at θ = 0). 


Figure <ref> presents results from   and illustrates how thermal conduction affects the evolution. It shows snapshots at t=2 of an adiabatic mountain (F = 0, left panel) and one evolved with a nonzero thermal flux (F≠ 0, right panel). The mountain grows from ≈ 1.3 × 10^3 cm to ≈ 3.0 × 10^3 cm. It is ≈ 4% taller for the run with conduction, but its density is lower (by a factor ≲ 5 near θ = 0.1 where an underdense column forms; see Sec. 3.2) everywhere except at θ≈ 0 and θ≈ 0.2. This is a consequence of the continuity equation (<ref>), which demands that an increase in height is met with an overall decrease in mass density. Aside from the underdense column, thermal conduction has the effect of driving matter towards the pole, where the density attains a maximum of ρ_max = 1.3 × 10^13 g cm^-3 for the run without conduction (left panel) and ρ_max = 1.5 × 10^13 g cm^-3 for the run with conduction (right panel) (≈ 16 % increase).   An analogous thermal softening phenomenon occurs in crustquake models, where thermal transport amplifies shear stresses felt in the neutron star crust <cit.>. 





Evolution with F≠ 0 tends to widen the magnetic field contours (cf. Fig. <ref>), because the mountain spreads and drags the field-lines with it [the time-dependent version of the flux-freezing condition (<ref>)]. Thermal conduction causes matter to be shifted both towards the pole and towards the base of the mountain, causing magnetic `pockets' to form near θ = 0 (see also Sec. 3.2 and Fig. 6, where they are clearer), as the field lines bend around the drifting matter. Overall, the magnetic field is weakened, going from a maximum strength of ||_max = 2.9 × 10^15 G to ||_max = 1.9 × 10^15 G and ||_max = 2.0 × 10^15 G at t= 2 without and with conduction, respectively. Away from the pole, the locations of the maxima and minima of |B| are largely unaffected by conduction.



















 §.§ Mass density and magnetic field evolution






In Figure <ref> we plot ρ contours and magnetic field lines for model A (top row) with M_a =  2.4 × 10^-5 M_⊙≈ 1.0 M_c, model B (middle row) with M_a = 3.0 × 10^-8 M_⊙≈ 1.2 M_c, and model C (bottom row) with M_a = 2.0 × 10^-6 M_⊙≈ 1.0 M_c for times t=0 (left panel) and t = 2 without (middle panel) and with (right panel) thermal conduction. The initial state is read from the Grad-Shafranov output and is the same for both runs for any given EOS.  All three runs have similar[Ideally we would keep M_a/M_c exactly the same across all runs. However, this is impractical because one does not know what μ is prior to running the Grad-Shafranov code for a given M_a, and it is computationally expensive to try and tune M_a/M_c exactly.] values of M_a/M_c∼ 1, implying that μ at t=0 is approximately halved by burial in all three models (see Sec. 4.1). We emphasise that models B and C are poor approximations to a realistic crust. They are included throughout mainly to give the reader a general sense of how the mountain structure depends on the EOS as well as to make contact with previous work for completeness.





The model A mountain (top row of Fig. <ref>) grows taller over time, thereby reducing ρ̅; we find = 3.3 × 10^3 cm for both F = 0 and F≠ 0 at t= 2, compared to = 1.2 × 10^3 cm at t=0. The magnetic pole remains the densest region with maximum densities ρ_max = 2.38 × 10^13 g cm^-3 at t=0 and ρ_max = 1.43 × 10^13 g cm^-3 at t= 2 without thermal conduction and ρ_max = 1.45 × 10^13 g cm^-3 at t= 2 with thermal conduction. A narrow, underdense (ρ≲ 4 × 10^12 g cm^-3) column forms at θ≈ 0.1 for the run with thermal conduction. Field lines are noticeably distorted from their t=0 state by flux freezing through a combination of poleward flow, which pushes them towards θ = 0, and stretching caused by the ∼ three-fold increase in , which pushes them radially outward.

Equation (<ref>) implies T ∝ρ^Γ - 1. Hence we expect thermal conduction to be less influential in model A (Γ = 1.18) than in models B and C (Γ = 5/3). The evolution of adiabatic model B (middle row of Fig. <ref>) is noticeably affected by the non-zero flux term F (cf. middle and right panels). The matter column near the pole has height (θ≈ 0) = 1.8 × 10^4 cm for F = 0 and (θ≈ 0) = 3.0 × 10^4 cm for F≠ 0. Matter concentrates more at the base of the mountain (r ≈ R_in) for F≠ 0, reaching peak densities of ρ_max = 6.8 × 10^8 g cm^-3 for F = 0 and ρ_max = 9.2 × 10^8 g cm^-3 for F≠ 0 at θ = 0. Both evolved mountains are denser than the initial state (ρ_max = 6.4 × 10^8 g cm^-3). Thermal conduction drives matter towards the pole, like what is seen in Fig. <ref>. The mountain grows taller, albeit comparatively less so than for model A, going from peak altitudes ≈ 4.0 × 10^4 cm to ≈ 4.2 × 10^4 cm with F = 0 and ≈ 4.3 × 10^4 cm for F≠ 0.







Adiabatic model C (bottom row of Fig. <ref>) evolves like model B. The mountain grows taller on the conduction time-scale, going from ≈ 4.0 × 10^3 cm at t=0 to ≈ 5.0 × 10^3 cm for F = 0 and ≈ 5.3 × 10^3 cm for F≠ 0. At t=0, the density maximum ρ_max = 3.7 × 10^11 g cm^-3 lies at r ≈ R_in and θ=0. After evolution the density reaches maximum values at the same location of ρ_max = 3.5 × 10^11 g cm^-3 (≈ 5% decrease) for F = 0 and ρ_max = 4.0 × 10^11 g cm^-3 (≈ 8% increase) for F≠ 0  at t=2. The compression of matter at the pole suggests that F≠ 0 acts to `soften' the effective EOS.











By inspecting the magnetic field lines in Fig. <ref>, we see that |B| evolves similarly to ρ due to flux freezing. As discussed above, thermal conduction drives matter towards the pole, shifting B accordingly <cit.>. Hence || decreases on the whole as time passes, most dramatically in the case of model C, which predicts ||_max = 9.7 × 10^14 G at t=0, ||_max = 6.5 × 10^14 G for F = 0 at t=2, and ||_max = 6.6 × 10^14 G for F≠ 0. This is similar to what occurs for the representative example discussed in Sec. 3.1 and tests with different grid resolutions (not plotted), suggesting the possibility that |B|_max increases slightly for runs with conduction, independent of the EOS. Additional, higher-resolution convergence tests (cf. Appendix A) can be undertaken, if future observational applications warrant. The formation of dense filamentary regions (ρ∼ 10^-1ρ_max) for runs with F≠ 0 at θ≈ 1.0 causes several magnetic `pockets' to form near the pole, as the magnetic field lines bend around poleward-drifting matter. The formation of filaments near the magnetic pole, as observed across all simulations with thermal conduction, may stem from thermal Parker-like instabilities <cit.>. These instabilities introduce `finger-like' density structures, which emerge due to the propagation of contact discontinuities between lighter and denser sections of fluid <cit.>. Although we only have one fluid in our model, the strong dependence of the conduction coefficient (<ref>) on the local magnetic field strength, which varies strongly near the pole, may cause this `fenced-off' behaviour. It has been shown that unstable modes grow faster in the presence of anisotropic thermal conduction <cit.>.













 §.§ Heat flux





Figure <ref> displays contours of temperature T (colour scale) and the direction of the thermal flux F (arrows) extracted from the runs performed in Figure <ref>. We seek to identify the existence, and evolution, of thermal hot spots [e.g. <cit.>]. 


The temperature varies gradually with r and θ for model A (top row), because the polytropic index Γ is nearly unity. Nevertheless, the maximum (at the pole) and minimum (at the mountain-atmosphere interface) values of T are in the ratio T_max/T_min≈ 10, a significant contrast. The thermal flux is predominantly directed towards the base of the mountain at all times, but becomes more `noisy' at large t, when local hot spots form. At t=2 we see that the temperature profile becomes more uniform away from the pole, suggesting that the model evolves towards an isothermal end state [d/dt(p/ρ) ≈ 0 for t ≫], even when thermal conduction is not implemented.  At θ≈ 0.1 for the run with F≠ 0 (right panel) we see a region of relatively low temperature (T_max/T ∼ 6) form. This `heat sink' is underdense as seen in Fig. <ref>, and is surrounded by the local hot spots (T_spots∼ 5 × 10^9 K) described above. 

For Model B (middle row) we have |F| ≈ 10^23 erg cm^-2s^-1, and F is predominantly directed towards the pole at θ =0. In the run without conduction (middle panel) we see that F is almost indistinguishable from its t=0 counterpart except near the equator where heat flows into a hot column (T ≈ 6 × 10^9 K) at θ = 1.3. The temperature evolves like the density, i.e. growing and spreading with ρ (cf. Fig. <ref>). When conduction is switched on, a hot region (T ≳ 8 × 10^9 K) forms near the pole which extends to the mountain-atmosphere interface at ≈ 4.2 × 10^4 cm. The flux is highest at the mountain-atmosphere interface and at altitude r-R_in≈ 1.3 × 10^4 cm. The flux is directed in different directions throughout the column, suggesting that localised hot spots form in the densest part of the mountain. Away from the pole (θ≳ 0.1), the heat flow is small (|F| ≲ 10^-2 |F|_max).





In Model C (right panel) at t=0, we see that heat flows towards the pole (|F| ≈ 10^21 erg cm^-2s^-1) and away from the equator (|F| ≈ 10^20 erg cm^-2s^-1). At t=2, however, heat flows from the top of the mountain to the base near the equator (θ≈ 1.3), and little heat (|F| ≲ 10^-4 |F|_max) flows near the pole, where a hot column (T ≈ 8 × 10^8 K) develops in a manner to similar to model B.  Because of relation (<ref>), the temperature profile evolves like the density and increases as the mountain grows (T ∝ρ^Γ-1); cf. Fig <ref>. Overall, the initially polytropic mountains respond similarly to thermal conduction by forming hot spots near the equator at θ≈ 1.3 (models B and C) and near the pole at θ≈ 0.1 (all models), where |F| is largest.







§ GLOBAL OBSERVABLES


In this section we consider the evolution on the conduction time-scale ≫ (<ref>) of the global observables μ (Sec. 4.1) and ϵ (Sec. 4.2) derived from runs of models A, B, and C. Simulations of mountains on the Alfvén time-scale  have been performed previously using the codes  <cit.> and  <cit.>.

 



 §.§ Magnetic dipole moment


The theory of magnetic burial predicts that the global magnetic dipole moment for an axisymmetric mountain,

    μ = 3 r^4/4∫^1_-1 d ( cosθ) cosθ B_r(r,θ),

evaluated at r ≳ R_m, decreases as a function of M_a <cit.>. In order to explore the relationship between burial, accreted mass, and thermal conduction for different EOS, we calculate μ from for a variety of  simulations with thermal conduction switched on.




Figure <ref> shows how μ (normalized to the pre-accretion value μ_i) evolves due to thermal conduction for various accreted masses (0.1 ≤ M_a/M_c≤ 1; left to right) and initially adiabatic EOS (1.18 ≤Γ≤ 5/3; top to bottom). Each panel displays μ(t) on a logarithmic temporal scale to capture both the MHD (t ≪) and thermal (t ≳) dynamics. Again we emphasise that model A (top row of Fig. <ref>) corresponds most closely to an astrophysically realistic accreted crust. Models B and C are included for completeness to illustrate EOS-related trends and make contact with previous work <cit.>.




All the mountains depicted in Fig. <ref> undergo an initially violent phase within ≲ 10^2, during which μ drops then rises. The behaviour observed in model A is similar to what is seen in Figures 6 and 14 of <cit.> for example. It is largely driven by the MHD reconfiguration of the mountain rather than thermal conduction (10^2≪)  <cit.>. We find that μ decreases slightly for all mountains (maximum of ≈ 7% for model A with M_a/M_c≈ 0.2), independent of the EOS, from t=0 to t ∼ 10^2, consistent with previous  simulations <cit.>. Note that the Grad-Shafranov equilibria, and hence the evolution, are insensitive to the exact value of the initial dipole moment μ_i provided that we have M_a / M_c≲ 10 [cf. the scaling law (1) introduced by <cit.>] <cit.>. In this context, insensitive means that μ/μ_i depends on μ_i only through the ratio M_a/M_c and not on μ_i in isolation. Since <cit.> found that M_c∝μ_i^2, the insensitivity condition M_a / M_c≲ 10 translates into an EOS-dependent lower bound for μ_i. For the astrophysically relevant model A, we require [see expression (B26) of <cit.>]

    ( μ_i/3.2 × 10^30 G cm^3)^2 ≳ 0.4 ( M_a/10^-4 M_⊙) ( R_⋆/10^6 cm)^6,

which is safely applicable to many LMXB systems, at least within the early stages of accretion <cit.>. For M_a / M_c≳ 10, the Grad-Shafranov modelling breaks down, and it is an open question whether the results are sensitive to μ_i or not.

On the longer time-scale t ≳, the behaviour of μ is qualitatively similar for all three initially polytropic EOS. In all cases, μ increases beyond μ_i; runs with M_a/M_c≳ 0.1 lead to μ(t ≳) > μ(t=0). In effect thermal conduction resurrects some of the buried field, e.g. μ increases by ≈ 14% in the case of model C with M_a/M_c≈ 0.2. For runs with M_a/M_c≳ 0.5, μ increases significantly from its initial value at t=0, e.g. by up to ≈ 80% in the case of model B for M_a = 3.0 × 10^-8 M_⊙≈ 1.2 M_c.

For a realistic accreted crust (model A), we find 0.64 ≤μ/μ_i≤ 0.83 for 0.2 ≤ M_a/M_c≤ 1 at t = 2. Small changes in μ for t ≳ in model A suggest that conduction plays a comparatively minor role in the evolution of astrophysically realistic mountains. Nevertheless we find that substantial (μ/μ_i≲ 0.5) magnetic burial requires significantly greater accreted masses than previously estimated by <cit.>; for example, μ(t=2)/μ(t=0) = 1.28 for M_a/M_c(t=0) ≈ 1 suggests an increase in the characteristic mass M_c at t = 2 by a factor ∼ 2. 



The inclusion of thermal conduction has the effect of partially resurrecting the buried field by increasing μ, which is similar to `softening' the EOS [as found by <cit.>]. The comparatively small increase in μ for model A (see Fig. <ref>) implies that the realistic EOS softens less than for models B and C. This is expected because the polytropic index Γ  =1.18 is closer to unity (i.e. nearly isothermal), implying that T is smaller than for the isentropic gas models B and C.













 §.§ Mass ellipticity

The characteristic gravitational wave strain emitted by a continuous-wave source is <cit.>

    h_c = (128 π^4/15)^1/2G I_zzν^2 |ϵ|/d c^2,



where I_jk is the moment-of-inertia tensor, ν is the spin frequency, d is the distance from the Earth to the source, and ϵ is the mass ellipticity,

    ϵ = I_xx - I_yy/I_zz.

The magnitude of ϵ represents the primary uncertainty in estimating h_c in practical astrophysics applications [see e.g. <cit.>; though cf. <cit.>]. Here we can calculate ϵ directly from (<ref>) using ρ as output by . Thus we can explore the effects of thermal conduction on the detectability of magnetic mountains using ground-based interferometers such as the Laser Interferometer Gravitational-Wave Observatory (LIGO) <cit.>.


Figure <ref> plots ϵ against time (in units of ) for different EOS and values of M_a. The layout is the same as in Fig. <ref>. All of the runs yield ϵ < 0, indicating that the star is prolate[Since <cit.> found that three-dimensional simulations of magnetic mountains relax to an almost axisymmetric state after a few Alfvén times (see Footnote 1), we expect the star to be prolate even without the assumption of axial symmetry.]; the mountain is densest at the magnetic pole <cit.>. 



In contrast to the magnetic dipole moment (Sec. 4.1), |ϵ| increases with time for t ≳ 10^2 for all runs with M_a/M_c≳ 0.1, by up to ≈ 45 % in the astrophysically realistic model A with M_a/M_c≈ 1. We also find smaller but still significant increases in |ϵ| in models B (≈ 30 % with M_a/M_c≈ 1), and C (≈ 27% with M_a/M_c≈ 1). This result is consistent with the leading-order behaviour of ϵ given by (<ref>), which implies that |ϵ| increases with M_c <cit.>. As noted in Sec. 4.1, all runs display a ≳ two-fold increase in M_c at t = 2. Thermal conduction tends to facilitate the poleward drift of matter (see Sec. 3.2), thereby making the star more prolate. 





The wobble angle of a precessing prolate star tends to grow, until the rotation and principal axes are orthogonal <cit.>, which is the optimal state for gravitational wave emission. Hence, an accreting neutron star with a prolate magnetic mountain may be harder to detect than an isolated magnetar with the same |ϵ|, which is oblate <cit.>. Note that, as discussed in Sec. 2.3, the values of |ϵ| presented in this section should be treated as upper limits since we do not model sinking <cit.>.

A summary of simulation parameters and results is given in Table <ref>. 








§ LONG-TERM THERMAL RELAXATION


In this section we explore the long-term thermal relaxation of a representative example of an astrophysically realistic mountain, namely model A with M_a = 1.8 × 10^-5 M_⊙≈ 0.58 M_c. Long-term (t ≫) simulations face numerical difficulties because of the wide range of time-scales in the problem. For a typical mountain, maintaining a resolution of 128 × 128 grid points (see Appendix A) requires a time-step satisfying Δ t ≲ 10^-7 to avoid numerical instabilities. It is impractical to evolve the simulation for long times ( t ≫). Lower-resolution runs (e.g. 64 × 64) fail catastrophically at t ≳, because steep gradients are handled poorly at the now `blurry' mountain-atmosphere interface; one ends up with ρ < 0 in places, for example. To circumvent these difficulties, we artificially increase the conduction coefficients κ_⊥ and κ_|| to accelerate thermal relaxation; <cit.> took a similar approach to accelerate Ohmic decay. Increasing κ by a factor ≳ 50 causes the super-time-stepping algorithm to fail, when the parabolic Courant condition is eventually violated [see Appendix A and <cit.>]. However, for an acceleration factor of ≲ 50, the simulation is stable. 

In order to increase the conductivities artificially, we set κ_⊥↦ξκ_⊥ and κ_||↦ξκ_||, where 1 ≤ξ≤ 40 is a dimensionless constant. Equations (<ref>) and (<ref>) imply ≈ 112/ξ s≈ 5.6 × 10^6/ξ. Figures <ref> and <ref> display μ and ϵ respectively as functions of time for ξ = 1,10,20,30,40. To read the horizontal-axis for the ξ = 30 case, for example, a value on the axis of 10^-1 implies that an interval lasting t=  3 has effectively elapsed. The longest run effectively extends over the interval 0 ≤ t ≤ 80. 



Both μ and |ϵ| increase with ξ. In other words, as κ_⊥ and κ_|| increase, magnetic burial is mitigated, while the gravitational wave strain increases. Thermal conduction pushes matter towards the pole (as in Sec. 3.2 and Fig. <ref>). Increasing F by a factor ξ amplifies polarward transport, i.e. ρ increases at the pole, as ξ increases, which is why |ϵ| increases with ξ and the star becomes more prolate. Increasing ρ near the pole effectively reduces the fraction p / ρ^Γ = k_Γ there; i.e. increasing ξ can be thought of as reducing the effective polytropic constant and `softening' the EOS by a factor related to ξ; cf. Table <ref>. Hence, initially adiabatic mountains evolved with high F come to resemble isothermal mountains at t ≫ [compare Fig. <ref> with Figure 8 of <cit.>]. For example, for ξ = 40, we have μ(t ≫)/μ_i = 0.73, cf. μ(t=0)/μ_i = 0.63. By comparison, we find μ/μ_i= 0.73 and μ/ μ_i = 0.63 from isothermal Grad-Shafranov simulations (softest EOS) for M_a≈ 10^-5 M_⊙ and M_a≈ 10^-4 M_⊙ respectively. Comparing the ξ = 1 and ξ = 40 final states, we find that μ and ϵ differ by ≤ 3% and ≤ 1% respectively. The trends discussed above are also evident in simulations with different EOS and grid resolutions; see Appendix A.



 


We see that μ and ϵ continue to oscillate after t ≳ 2 albeit with small amplitude (≲ 1% peak to peak). The fluctuations at t ≫ persist, because some parts of the mountain take longer to settle down than others. In reality, heat transport occurs more slowly than average in cold regions (T_e≪ 10^8 K), meaning that conduction continues to affect cold parts of the mountain (whose effective conduction time-scales are longer than the volume-averaged value ), even after the rest of the mountain relaxes thermally. These cold regions, however, do not play a dominant role in determining μ or ϵ, as the density is low there. 












§ CONCLUSIONS




In this paper we explore the effects of thermal conduction on the evolution of accretion-built mountains on neutron stars for time-scales t ≳ (Secs. 3 and 4) and t ≫ (Sec. 5) using the MHD code  <cit.>. The initial states are generated from the Grad-Shafranov equation for a range of initially polytropic EOS documented in Table <ref>  <cit.>. Model A approximates a realistic, non-catalysed, accreted crust for densities in the range 10^8≤ρ /  g cm^-3≤ 10^14 <cit.>. Models B and C approach the realistic EOS in the low-ρ and neutron drip density regimes respectively, and are included for completeness to illustrate EOS-related trends and make contact with previous work. 
The theory of magnetic burial predicts that, as matter piles up on the stellar surface, the dipole moment is reduced in accreting neutron stars in accord with the observed μ versus M_a relations, e.g. <cit.>. We find that thermal conduction has the effect of pushing accreted matter back towards the magnetic pole, where T is greatest, thereby partially resurrecting the buried field and increasing μ while making the star more prolate. On the conduction time-scale, we find a quasi-static increase in the mountain's characteristic mass M_c [defined above equation (<ref>)] starting from an adiabatic initial state. 
Hence achieving a given μ/μ_i value requires higher M_a, in general, than estimated by <cit.>. The main trends are summarised in Table <ref>. 


Gravitational radiation back-reaction can stall the spinup of the neutron star in a low-mass X-ray binary (LMXB) at hectohertz frequencies <cit.>, explaining the observation that LMXBs spin slower (ν_s≲ 650  Hz) than otherwise expected <cit.>. The results in this paper suggest that the effective EOS of mountain matter may be softer than previously estimated, when thermal conductivity is included, leading to a proportionally higher gravitational wave strain (<ref>). This strengthens the argument for targeting LMXBs such as Sco X-1 for searches with facilities like LIGO <cit.>. The increase in |ϵ| combined with the decrease in μ/μ_i at t ≳ for all runs with M_a/M_c≳ 0.1 performed in this paper suggests that stars with significantly buried (μ≪μ_i) magnetic fields may prove better gravitational wave candidates than previous estimates indicated <cit.>. However, we stress that the systematic and numerical (see Appendix A) uncertainties present within our models suggest that the effects of thermal conduction are likely to be small compared to other physical effects not implemented here, such as sinking <cit.>. 


In addition to searching for gravitational waves and measuring the global dipole moment, one can test the magnetic burial scenario by studying type I X-ray bursts <cit.>. For adiabatic initial states, we find that mountains develop hot (T ≳ 10^9 K) spots with large thermal fluxes near both the pole (θ = 0.1) and the equator (θ = 1.3) for a wide range of accreted masses (see Table <ref> and Fig. <ref>). Dense filamentary regions also develop, especially for M_a/M_c≳ 1. These effects cooperate to produce localized hot patches “fenced off" by intense magnetic fields, whose number increases with M_a [cf. <cit.>]. The hot spots may individually provide fuel for type I X-ray bursts which do not spread across the entire stellar surface, if the magnetic fences are intense enough to inhibit cross-field thermal transport <cit.>. X-ray observations of significant heat fluxes near the magnetic pole of a neutron star in an LMXB, as broadly predicted by our simulations, may be related to the magnetic mountain physics <cit.>. Thermal fluxes out of the hot spots may also amplify shear stresses felt by the neutron star crust <cit.>. A detailed analysis of hot-spot phenomena and their observational consequences will be conducted in future work. Another avenue to probe accretion mound physics comes from cyclotron features <cit.>. <cit.> showed that one can discriminate, in principle, between magnetic mountain properties (e.g. EOS) by studying the line energy, width, and depth of theoretical cyclotron resonant scattering features from accreting neutron stars. These cyclotron features are, however, unlikely to be detected in the near future as it requires further development of sensitive X-ray polarimeters <cit.>. 



















§ ACKNOWLEDGMENTS

We thank Maxim Priymak and Donald Payne for permission to modify and use the Grad-Shafranov solver. We thank Dipanjan Mukherjee for expert instruction on the use of . We thank HoChan Cheon for designing an early version of the script that converts Grad-Shafranov output into a format suitable for PLUTO input. We thank Patrick Clearwater, Brynmor Haskell, and Alpha Mastrano for discussions. We thank the anonymous referee for their carefully considered comments, and for providing the Skyrme EOS data used in Sec. 2.3 to estimate sinking depths. This work was supported in part by an Australian Postgraduate Award, the Albert Shimmins fund, and an Australian Research Council Discovery Project grant. 






[Aasi et al.(2014)]aasi14 Aasi, J., Abadie, J., Abbott, B. P., et al. 2014, , 785, 119 

[Abbott et al.(2007)]scox1 Abbott, B., Abbott, R., Adhikari, R., et al. 2007, , 76, 082001 



[Akgün et al.(2018)]akgun18 Akgün, T., Cerdá-Durán, P., Miralles, J. A., & Pons, J. A. 2018, , 481, 5331 

[Alexiades et al.(1996)]aag96 Alexiades, V., Amiez, G., & Gremaud, P. 1996, Com. Num. Meth. Eng., 12, 31


[Alfvén(1943)]alf43 Alfvén, H. 1943, Arkiv for Astronomi, 29, 1 





[Arons(1993)]arons93 Arons, J. 1993, , 408, 160 

[Arzoumanian et al.(2002)]arzo Arzoumanian, Z., Chernoff, D. F., & Cordes, J. M. 2002, , 568, 289 

[Balbus(1986)]bal86 Balbus, S. A. 1986, , 304, 787 

[Balescu(1988)]balescu1 Balescu, R. 1988, Transport Processes in a Plasma (Amsterdam: North-Holland)

[Becker & Truemper(1997)]thermfluxpaper Becker, W., & Truemper, J. 1997, , 326, 682 

[Beckers(1992)]bec92 Beckers, J. M. 1992, SIAM Journal on Numerical Analysis, 29, 701

[Beloborodov & Levin(2014)]belo14 Beloborodov, A. M., & Levin, Y. 2014, , 794, L24 

[Bhattacharyya & Strohmayer(2006)]bhat06 Bhattacharyya, S., & Strohmayer, T. E. 2006, , 641, L53 

[Bildsten(1998)]bildstenpaper Bildsten, L. 1998, , 501, L89 

[Bisnovatyĭ-Kogan & Chechetkin(1979)]kogan79 Bisnovatyĭ-Kogan, G. S., & Chechetkin, V. M. 1979, Soviet Physics Uspekhi, 22, 89 

[Blondin & Freese(1986)]blon86 Blondin, J. M., & Freese, K. 1986, , 323, 786 

[Brady et al.(1998)]brady Brady, P. R., Creighton, T., Cutler, C., & Schutz, B. F. 1998, , 57, 2101 

[Braginskii(1965)]braginskii65 Braginskii, S. I. 1965, Reviews of Plasma Physics, 1, 205 

[Brown(2000)]brown00 Brown, E. F. 2000, , 531, 988 

[Brown & Bildsten(1998)]bil98 Brown, E. F., & Bildsten, L. 1998, , 496, 915 

[Brown et al.(1998)]rutledge Brown, E. F., Bildsten, L., & Rutledge, R. E. 1998, , 504, L95 


[Čada & Torrilhon(2009)]lim03 Čada, M., & Torrilhon, M. 2009, Journal of Computational Physics, 228, 4118 

[Cavecchi et al.(2017)]watts17 Cavecchi, Y., Watts, A. L., & Galloway, D. K. 2017, , 851, 1 

[Chakrabarty et al.(2003)]chaknat03 Chakrabarty, D., Morgan, E. H., Muno, M. P., et al. 2003, , 424, 42 

[Chamel & Haensel(2008)]cham08 Chamel, N., & Haensel, P. 2008, Living Reviews in Relativity, 11, 10 

[Chamel et al.(2015)]cham15 Chamel, N., Fantina, A. F., Zdunik, J. L., & Haensel, P. 2015, , 91, 055803 

[Chandrasekhar(1956)]C56 Chandrasekhar, S. 1956, , 124, 232 

[Chandrasekhar(1967)]ch76 Chandrasekhar, S. 1967, New York: Dover, 1967,  

[Chugunov & Horowitz(2010)]chug10 Chugunov, A. I., & Horowitz, C. J. 2010, , 407, L54 

[Courant & Hilbert(1953)]lagchap Courant, R., & Hilbert, D. 1953, Methods of Mathematical Physics, Vol. I (New York: Interscience Publication)

[Cumming(2004)]cumm04 Cumming, A. 2004, IAU Colloq. 190: Magnetic Cataclysmic Variables, 315, 58 


[Cumming et al.(2001)]cumm01 Cumming, A., Zweibel, E., & Bildsten, L. 2001, , 557, 958 


[Cutler(2002)]cutler02 Cutler, C. 2002, , 66, 084025 

[Douchin & Haensel(2001)]douchin01 Douchin, F., & Haensel, P. 2001, , 380, 151 

[Faucher-Giguère & Kaspi(2006)]kaspi Faucher-Giguère, C. A., & Kaspi, V. M. 2006, , 643, 332 

[Field(1965)]field1965 Field, G. B. 1965, , 142, 531 

[Fujimoto et al.(1984)]fuji84 Fujimoto, M. Y., Hanawa, T., Iben, I., Jr., & Richardson, M. B. 1984, , 278, 813 

[Galloway et al.(2008)]gal08 Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, , 179, 360-422 

[Galloway et al.(2014)]gal14 Galloway, D. K., Premachandra, S., Steeghs, D., et al. 2014, , 781, 14 

[Gardiner & Stone(2005)]godunov05 Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509 

[Geppert & Viganò(2014)]gv14 Geppert, U., & Viganò, D. 2014, , 444, 3198

[Glampedakis & Andersson(2006)]glam06 Glampedakis, K., & Andersson, N. 2006, , 74, 044040 

[Goldreich & Reisenegger(1992)]reissen1 Goldreich, P., & Reisenegger, A. 1992, , 395, 250 



[Harry & LIGO Scientific Collaboration(2010)]abb08 Harry, G. M., & LIGO Scientific Collaboration 2010, Classical and Quantum Gravity, 27, 084006 

[Haskell et al.(2006)]haskell06 Haskell, B., Jones, D. I., & Andersson, N. 2006, , 373, 1423 

[Haskell et al.(2015)]hask15 Haskell, B., Priymak, M., Patruno, A., et al. 2015, , 450, 2393 

[Haensel & Zdunik(1990a)]haen89 Haensel, P., & Zdunik, J. L. 1990, , 227, 431 

[Haensel & Zdunik(1990b)]haen90 Haensel, P., & Zdunik, J. L. 1990, , 229, 117 


[van den Heuvel & Bitzaraki(1995)]hue95 van den Heuvel, E. P. J., & Bitzaraki, O. 1995, , 297, L41 

[Johnson-McDaniel & Owen(2013)]mcdaniel13 Johnson-McDaniel, N. K., & Owen, B. J. 2013, , 88, 044004 

[Keek et al.(2010)]gal10 Keek, L., Galloway, D. K., in't Zand, J. J. M., & Heger, A. 2010, , 718, 292 

[Klimchuk & Sturrock(1989)]sturrock89 Klimchuk, J. A., & Sturrock, P. A. 1989, , 345, 1034 


[Konar & Choudhuri(2002)]kon02 Konar, S., & Choudhuri, A. R. 2002, Bulletin of the Astronomical Society of India, 30, 697 

[Konar & Choudhuri(2004)]kon04 Konar, S., & Choudhuri, A. R. 2004, , 348, 661 

[Konar(2010)]kon11 Konar, S. 2010, , 409, 259 

[Kosiński & Hanasz(2006)]kos06 Kosiński, R., & Hanasz, M. 2006, Astronomische Nachrichten, 327, 479 


[Kundu & Cohen(2008)]kundu Kundu, P. K., & Cohen, I. M. 2008, Fluid Mechanics: Fourth Edition (London, England: Academic Press)

[Landau & Lifshitz(1959)]LL59 Landau, L. D., & Lifshitz, E. M. 1959, Fluid Mechanics (Oxford: Pergamon Press)

[Lasky(2015)]lasky15 Lasky, P. D. 2015, PASA, 32, e034 



[Lecoanet et al.(2012)]leco12 Lecoanet, D., Parrish, I. J., & Quataert, E. 2012, , 423, 1866 

[Litwin et al.(2001)]litwin01 Litwin, C., Brown, E. F., & Rosner, R. 2001, , 553, 788 

[Mackie & Baym(1977)]mackie77 Mackie, F. D., & Baym, G. 1977, Nuclear Physics A, 285, 332 

[Mastrano et al.(2011)]mmra11 Mastrano, A., Melatos, A., Reisenegger, A., & Akgün, T. 2011, , 417, 2288 

[Mastrano & Melatos(2012)]mm12 Mastrano, A., & Melatos, A. 2012, , 421, 760 

[Mastrano et al.(2015)]MSM15 Mastrano, A., Suvorov, A. G., & Melatos, A. 2015, , 447, 3475 

[Melatos & Payne(2005)]mp05 Melatos, A., & Payne, D. J. B. 2005, , 623, 1044 

[Melatos & Phinney(2001)]melphi01 Melatos, A., & Phinney, E. S. 2001, PASA, 18, 421 

[Mignone et al.(2007)]pluto Mignone, A., Bodo, G., Massaglia, S., et al. 2007, , 170, 228 

[Mignone & Tzeferacos(2010)]hdivcl Mignone, A., & Tzeferacos, P. 2010, Journal of Computational Physics, 229, 2117 

[Miralda-Escude et al.(1990)]miralda90 Miralda-Escude, J., Paczynski, B., & Haensel, P. 1990, , 362, 572 


[Misanovic et al.(2010)]misanovic10 Misanovic, Z., Galloway, D. K., & Cooper, R. L. 2010, , 718, 947 

[Miyoshi & Kusano(2005)]MK05 Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315 

[Mouschovias(1974)]mous74 Mouschovias, T. C. 1974, , 192, 37 

[Mouschovias et al.(2009)]mous09 Mouschovias, T. C., Kunz, M. W., & Christie, D. A. 2009, , 397, 14 

[Mukherjee & Bhattacharya(2012)]dip12 Mukherjee, D., & Bhattacharya, D. 2012, , 420, 720 

[Mukherjee et al.(2013a)]dip13a Mukherjee, D., Bhattacharya, D., & Mignone, A. 2013, , 430, 1976 

[Mukherjee et al.(2013b)]dip13b Mukherjee, D., Bhattacharya, D., & Mignone, A. 2013, , 435, 718 



[Narayan & Heyl(2003)]heyl03 Narayan, R., & Heyl, J. S. 2003, , 599, 419 



[Nishimura(2005)]nish05 Nishimura, O. 2005, , 57, 769 

[Parker(1953)]park1953 Parker, E. N. 1953, , 117, 431 

[Patruno(2012)]pat12 Patruno, A. 2012, , 753, L12 

[Payne & Melatos(2004)]pm04 Payne, D. J. B., & Melatos, A. 2004, , 351, 569 

[Payne & Melatos(2006a)]pm06 Payne, D. J. B., & Melatos, A. 2006, , 641, 471 

[Payne & Melatos(2006b)]pmtherm Payne, D. J. B., & Melatos, A. 2006, , 652, 597 

[Payne & Melatos(2007)]pstab07 Payne, D. J. B., & Melatos, A. 2007, , 376, 609



[Potekhin(1999)]pot99 Potekhin, A. Y. 1999, , 351, 787 

[Potekhin et al.(2013)]pot13 Potekhin, A. Y., Fantina, A. F., Chamel, N., Pearson, J. M., & Goriely, S. 2013, , 560, A48 

[Press et al.(1986)]levmac Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Cambridge: University Press, 1986,  

[Priymak et al.(2014)]pri14 Priymak, M., Melatos, A., & Lasky, P. D. 2014, , 445, 2710 

[Priymak et al.(2011)]pri11 Priymak, M., Melatos, A., & Payne, D. J. B. 2011, , 417, 2696 

[Riles(2013)]riles13 Riles, K. 2013, Progress in Particle and Nuclear Physics, 68, 1 



[Sato(1979)]sato79 Sato, K. 1979, Progress of Theoretical Physics, 62, 957 

[Schatz et al.(1999)]schatz99 Schatz, H., Bildsten, L., Cumming, A., & Wiescher, M. 1999, , 524, 1014 

[Shapiro & Teukolsky(1983)]ST83 Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars: The physics of compact objects (New York: Wiley-Interscience)

[Shibazaki et al.(1989)]shibaz Shibazaki, N., Murakami, T., Shaham, J., & Nomoto, K. 1989, , 342, 656 



[Srinivasan et al.(1990)]srini90 Srinivasan, G., Bhattacharya, D., Muslimov, A. G., & Tsygan, A. J. 1990, Current Science, 59, 31 

[Stone & Gardiner(2007)]stone07 Stone, J. M., & Gardiner, T. 2007, , 671, 1726 

[Strohmayer & Bildsten(2003)]stroh03 Strohmayer, T., & Bildsten, L. 2006, Cambridge Astrophysics Series, 39, 113

[Suvorov et al.(2016a)]SMM16 Suvorov, A. G., Mastrano, A., & Melatos, A. 2016, , 456, 731 

[Suvorov et al.(2016b)]geppert16 Suvorov, A. G., Mastrano, A., & Geppert, U. 2016, , 459, 3407 

[Suvorov(2018)]suvorov18 Suvorov, A. G. 2018, , 98, 084026 

[Taam & van den Heuvel(1986)]taam86 Taam, R. E., & van den Heuvel, E. P. J. 1986, , 305, 235 

[Thorne(1980)]thorne80 Thorne, K. S. 1980, Reviews of Modern Physics, 52, 299 

[Urpin & Geppert(1995)]urp95 Urpin, V., & Geppert, U. 1995, , 275, 1117 

[Ushomirsky et al.(2000)]usho00 Ushomirsky, G., Cutler, C., & Bildsten, L. 2000, , 319, 902 

[Vigelius & Melatos(2008)]vigstab Vigelius, M., & Melatos, A. 2008, , 386, 1294 

[Vigelius & Melatos(2009a)]vig09 Vigelius, M., & Melatos, A. 2009, , 395, 1963 

[Vigelius & Melatos(2009b)]vigohm09 Vigelius, M., & Melatos, A. 2009, , 395, 1985 

[Vigelius & Melatos(2009c)]vig09c Vigelius, M., & Melatos, A. 2009, , 395, 1972 

[Wang et al.(2012)]wang12 Wang, J., Zhang, C. M., & Chang, H.-K. 2012, , 540, A100 

[Wette et al.(2010)]wette10 Wette, K., Vigelius, M., & Melatos, A. 2010, , 402, 1099 

[Yoshida(2013)]yoshida13 Yoshida, S. 2013, , 435, 893 

[Zdunik et al.(1992)]zdun92 Zdunik, J. L., Haensel, P., Paczynski, B., & Miralda-Escude, J. 1992, , 384, 129 

[Zhang & Kojima(2006)]zhang06 Zhang, C. M., & Kojima, Y. 2006, , 366, 137 









§ PLUTO SIMULATIONS


Complete documentation for the  code was published by <cit.>. The specific features we rely upon and optimize are discussed below.



 



 §.§ Grid and time step


We employ a static, two-dimensional, polar grid with  N_r× N_θ = 128 × 128 grid points. The radial grid comprises a logarithmic section with 100 points for R_in≤ r ≤ 2 R_mountain, and a uniform section with 28 points for 2 R_mountain≤ r ≤ R_m, where R_mountain is defined arbitrarily at t=0 as the innermost radial grid point with ρ≤ 10^-7ρ_max. A mixed grid captures features with sharply different length-scales and minimizes the interpolation errors discussed in Sec. 3. We find that including additional grid points in the atmosphere (R_mountain≤ r ≤ R_m) increases the computational cost without modifying perceptibly the observables computed in Sec. 4. The angular grid is uniformly spaced in θ over 0 ≤θ≤π/2.

We employ a Runge-Kutta third-order time-stepper for safety, although we find by experimentation that the results are essentially indistinguishable from the second-order variant. We employ the third-order finite-volume spatial integrator `Lim03' to interpolate between grid points <cit.>. This scheme resolves local minima with high precision, e.g. strong gradients at the mountain-atmosphere interface. We use a time-step Δ t = 5 × 10^-8 where  is determined through equation (<ref>). We print output files at various fractions of  depending on the specifics of the run; cf. the horizontal-axes on Figs. <ref> and <ref>. This Δ t is small enough to avoid Courant-Friedrichs-Lewy (CFL) instabilities for each run.

In order to avoid numerical instabilities we simulate the mountain together with an atmosphere which has a small but non-zero density taken as ρ_atm = 10^-8ρ_max. The atmosphere alleviates numerical difficulties associated with strong gradients and prevents the density from dipping below zero due to numerical fluctuations. We find that varying ρ_atm in the range 10^-9≤ρ_atm/ ρ_max≤ 10^-5 does not modify the observables discussed in Sec. 3 and 4 by more than ∼ 1%. The code crashes for ρ_atm≲ 10^-9ρ_max, because ρ dips below zero somewhere unless we set Δ t ≲ 10^-11, which is too expensive computationally. We define a flag that sets max(ρ,ρ_atm) ↦ρ at every grid point after each time step Δ t so that the atmosphere has a minimum density of ρ_atm for all t.

We find that the bilinear interpolation algorithm of  introduces errors of ≤ 1 % in the quantities calculated in Secs. 3 and 4 between the t=0 (grid-realigned)  output and the raw Grad-Shafranov data; cf. Figures 2 and 4 of <cit.> with Figs. <ref> and <ref>.




 §.§ Divergence cleaning and thermal conduction


Maxwell's equations require · = 0 at all times. Various strategies can be employed to minimise  numerical deviations from · = 0. For example, the extended hyperbolic divergence cleaning algorithm <cit.> introduces Lagrange multipliers into Faraday's law (<ref>). Inspection of  output files confirms that · vanishes to floating-point precision as a consequence of using this algorithm. The divergence cleaning algorithm is coupled with the approximate Riemann solver `hlld' designed to resolve shocks and strong gradient phenomena <cit.>.

Thermal conduction (see Sec. 3) is implemented via the super-time-stepping algorithm available in  and described in <cit.>. The energy equation (<ref>) has a parabolic Courant number C_p associated with it, which depends on the value of the conduction coefficients κ_|| and κ_⊥. Together with the usual Courant number condition <cit.>, we require C_p≤ 1/N_dim = 1/2 to avoid instabilities <cit.>. Super-time-stepping allows for flux terms to be treated in a separate `super-step' using operator splitting methods, so that Δ t need not be reduced to avoid parabolic CFL instabilities.




 §.§ Implementation, stability, and convergence tests


We test our   simulations in four ways. For implementation: (i) We compute the total mass of the simulation at each time-step to check for mass leakage. (ii) We check the surface dipole moment and the velocity field to ensure that the boundary conditions described in Sec. 2.2 are implemented faithfully. For convergence: (iii) We vary the grid parameters N_r and N_θ to check if the results depend on the spatial resolution (see Figs. <ref> and <ref>). For stability: (iv) We vary the super-time-stepping parameters, the CFL parameters, and the thermal conduction coefficients (i.e. checking if the conduction and no-conduction runs match smoothly in the limit κ→ 0).  The convergence of the Grad-Shafranov code described in Sec. 2.2 is studied fully by <cit.> and <cit.>. 

In Figures <ref> and <ref> we show the evolution of the ellipticity and dipole moment, respectively, for model B with M_a = 2 × 10^-8 M_⊙ and thermal conduction switched on, with N_r× N_θ = 96 × 96 grid points and varying values of ξ (this parameter is introduced to artificially scale the conduction coefficients, see Sec. 5). Two major points are evident from these plots. First, the trends associated with increasing ξ for model B are the same as was observed for model A in Sec. 5; |ϵ| and μ/μ_i are monotonically increasing with increasing values of ξ at t ≳τ_c, independent of the EOS and grid resolution. The second point concerns the convergence test (iii) detailed above: for the ξ = 1 run (black, diamonds), all simulation parameters are identical to those for the simulations performed in Sec. 4 for the same accreted mass and EOS (middle figure of the middle panel), except that the resolution is lower for the runs presented here. Comparing the final ellipticity and μ values from Figs. <ref> and <ref> with those presented in Table 3 for the higher resolution run, we see only a small (≲ 10 %) disparity at late times, with ϵ(96 × 96) = 1.65 × 10^-8, ϵ(128 × 128) = 1.44 × 10^-8, μ(96 × 96)/μ_i = 0.83, and μ(128 × 128)/μ_i = 0.87.







Figure <ref> plots the total mass as a function of time without (black, triangles) and with (blue, circles) thermal conduction for model B with M_a = 2 × 10^-8 M_⊙. We see that, after an initial adjustment phase, the total mass returns to M_a within ∼ 3 % (∼ 4 %) without (with) thermal conduction. This adjustment phase (t ≪τ_c) occurs for two separate reasons. The first is due to the artificial atmospheric density ρ_atm = 10^-8ρ_max, introduced to ensure that the simulation does not produce ρ < 0 at any point throughout the evolution. Some of this atmospheric mass actually gets pulled down into the mountain, after which the atmosphere resets, thus increasing the overall mass of the simulation slightly. Additionally, the Grad-Shafranov equilibria are defined over grids which are slightly different to those in PLUTO. Hence, at t=0, the MHD equations are not exactly satisfied in PLUTO, leading to a temporary increase in the total mass. These two effects combine to increase the total mass in the initial stages of evolution. Table A2 in <cit.> reports similar total mass changes during the adjustment phase. Figure <ref> is typical for runs performed in this paper.




In Figure <ref> we plot the surface dipole moment μ_S [equation (<ref>) evaluated at r = R_in] as a function of time without (black, triangles) and with (blue, circles) thermal conduction for model A with M_a = 1.8 × 10^-5 M_⊙. If the boundary conditions at the stellar surface [namely ψ(R_in,θ) = ψ_⋆sin^2θ] are implemented without numerical error, μ_S should keep its initial value μ_i. We see a slight variation (maximum of ∼ 2 %). Figure <ref> is typical for runs performed in this paper.





§ IDEAL-GAS APPROXIMATION TO THE EQUATION OF STATE


Strictly speaking, the accreted matter in the crust is partially degenerate <cit.>. In this appendix, we verify that it is reasonable to approximate the EOS by the ideal-gas formula (<ref>), for ease of use in , when calculating the perturbations to the mountain structure caused by thermal transport. The equilibrium configuration of the mountain before thermal transport is switched on is calculated for the full, degenerate, polytropic EOS (see Sec. 2).

In a Fermi-Dirac distribution, the mean occupancy n for a single-particle orbital with energy E is given by

    n(E) = 1/exp[E - σ/k_B T] +1 ,

where σ is the chemical potential, which is a function of ρ and T (in general). In the limit T → 0, n(E) tends to either 1 or 0 for E < σ or E > σ, respectively. The Fermi temperature T_F is defined through the chemical potential via

    σ = k_B T_F.

The dependence of σ on ρ and T is determined by integrating the mean occupancy to obtain the total particle number,

    N = ∫^∞_0 n(E) D(E) d E,

where D(E) is the density of states.

The pressure p is defined via the first law of thermodynamics, viz. 

    p =  -(T - ε S - σ N )/V,

where ε and S denote the internal energy and entropy, respectively, and -pV is the grand canonical potential. Substituting (<ref>) into the integral (<ref>) allows one to express T in terms of p and ρ for T ≠ 0, i.e. defines the EOS. One finds

    p m/ρ k_B T = ℱ(5/2,z)/ℱ(3/2,z),

for an arbitrary Fermi gas, with

    ℱ(ν,z) = 1/Γ(ν)∫^∞_0x^ν-1/z^-1 e^x +1  dx,

and fugacity z = exp( σ / k_B T) [see e.g. <cit.> for details]. Expression (<ref>) is plotted in Fig. <ref>.

In regions with T ≫ T_F, expression (<ref>) approaches the ideal gas law (<ref>),

    p ∝ρ k_B T.
 
In regions with T ≪ T_F, expression (<ref>) approaches a polytropic EOS <cit.>,

    p ∝ρ^5/3.

An accurate description of a realistic accreted crust lies between these two extremes <cit.>. The latter limit (<ref>) coincides with the Grad-Shafranov t=0 initial condition for degenerate, single-particle fluids, e.g. models B and C in Table <ref> and Figure <ref>. The Grad-Shafranov equilibria, calculated with (<ref>), adjust modestly, when thermal transport is switched on in   with (<ref>), suggesting that the equilibrium starting-point is broadly consistent with both (<ref>) and (<ref>), for the values of κ_|| and κ_⊥ relevant here.



As noted throughout the body of the paper, it is the hot regions of the mountain where thermal conduction modifies the hydrodynamic structure the most as time passes. This is expected because the (dominant) parallel thermal conductivity κ_|| scales as κ_||∝ T^5/2 through (<ref>). In hot regions, we have T ≳ T_F. Hence Fig. <ref> implies a ≲ 20% departure in p m/ ρ k_B T from the ideal gas law.
 
In Figure <ref> we plot contours of T_exact/T_ideal [i.e. T_exact from (<ref>) divided by T_ideal from (<ref>)] (left panel) and κ_||,exact / κ_||,ideal (similarly defined, right panel) for the realistic accreted crust model A with M_a = 1.8 × 10^-5 M_⊙≈ 0.58 M_c at t=0. We find T_exact/T_ideal≳ 0.85 throughout the bulk of the mountain, i.e. the temperature is overestimated by ≲ 15 % in the densest regions of the mountain, where most mass resides, for this representative simulation. This translates into a ≲ 40 % overestimate in κ_||≫κ_⊥ everywhere except at the mountain-atmosphere interface, where there is little mass, and the model breaks down anyway because of the artificial ρ_atm.





