Asymmetric accretion and thermal `mountains' in magnetized neutron star crusts

TL;DR

This paper investigates how magnetic field-induced asymmetries in accreting neutron star crusts can generate gravitational waves, with implications for star spin evolution and gravitational wave detection.

Contribution

It provides a realistic calculation of crustal asymmetries caused by magnetic fields and assesses their potential to produce observable gravitational waves.

Findings

Outer crust layers are strongly perturbed by magnetic fields.

Estimated quadrupoles are generally too small to explain observed spin-downs.

Strong shallow heating can enhance perturbations, potentially explaining observed phenomena.

Abstract

Accreting neutron stars are one of the main targets for continuous gravitational wave searches, as asymmetric accretion may lead to quadrupolar deformations, or `mountains', on the crust of the star, which source gravitational wave emission at twice the rotation frequency. The gravitational wave torque may also impact on the spin evolution of the star, possibly dictating the currently observed spin periods of neutron stars in Low Mass X-ray Binaries and leading to the increased spindown rate observed during accretion in PSR J1023+0038. Previous studies have shown that deformed reaction layers in the crust of the neutron star lead to thermal and compositional gradients that can lead to gravitational wave emission. However, there are no realistic constraints on the level of asymmetry that is expected. In this paper we consider a natural source of asymmetry, namely the magnetic field, and calculate the density and pressure perturbations that are expected in the crust of accreting neutron stars. In general we find that only the outermost reaction layers of the neutron star are strongly perturbed. The mass quadrupole that we estimate is generally small and cannot explain the increase of spin-down rate of PSR J1023+0038. However, if strong shallow heating sources are present at low densities in the crust, as cooling observations suggest, these layers will be strongly perturbed and the resulting quadrupole could explain the observed spindown of PSR J1023+0038, and lead to observable gravitational wave signals from systems with higher accretion rates.

Tables

Table 1. Table 1: Reactions in the crust of an accreting NS assuming that the ashes of X-ray bursts consist of pure Fe 56 superscript Fe 56 {}^{56}{\rm Fe} using the EDF BSk21 as described in Fantina et al. ( 2018 ) . We provide pressure P 𝑃 P , density ρ 𝜌 \rho , threshold energy E t ​ h subscript 𝐸 𝑡 ℎ E_{th} and heat release per accreted baryon Q M subscript 𝑄 𝑀 Q_{M} . For completeness we also provide the free neutron fraction X n subscript 𝑋 𝑛 X_{n} and the density jump in the layer Δ ​ ρ / ρ Δ 𝜌 𝜌 \Delta\rho/\rho .

$P$	$\rho$	Reactions	$X_n$	$\Delta \rho /\rho$	$E_{th}$	$Q_M$
(dyn cm-2)	(g cm-3)				(MeV)	(keV)
$6.48\times {10}^{26}$	$1.38\times {10}^9$	${}_{}{}^{56}\mathrm{Fe}\to {}_{}{}^{56}\mathrm{Cr}-2e^{-}+2{\nu }_e$	0	0.08	4.47	37.0
$1.83\times {10}^{28}$	$1.81\times {10}^{10}$	${}_{}{}^{56}\mathrm{Cr}\to {}_{}{}^{56}\mathrm{Ti}-2e^{-}+2{\nu }_e$	0	0.09	10.22	41.2
$1.06\times {10}^{29}$	$7.37\times {10}^{10}$	${}_{}{}^{56}\mathrm{Ti}\to {}_{}{}^{56}\mathrm{Ca}-2e^{-}+2{\nu }_e$	0	0.10	15.83	39.1
$3.43\times {10}^{29}$	$1.96\times {10}^{11}$	${}_{}{}^{56}\mathrm{Ca}{\to }^{56}\mathrm{Ar}-2e^{-}+2{\nu }_e$	0	0.11	21.22	8.1
$8.75\times {10}^{29}$	$4.38\times {10}^{11}$	${}_{}{}^{56}\mathrm{Ar}\to {}_{}{}^{55}\mathrm{Cl}+n-e^{-}+{\nu }_e$	0	0.06	26.55	0
$9.40\times {10}^{29}$	$4.79\times {10}^{11}$	${}_{}{}^{55}\mathrm{Cl}\to {}_{}{}^{53}S+\Delta N\cdot n-e^{-}+2{\nu }_e$	0.05	0.06	27.04	0
$1.18\times {10}^{30}$	$6.04\times {10}^{11}$	${}_{}{}^{53}S\to {}_{}{}^{48}\mathrm{Si}+\Delta N\cdot n-2e^{-}+2{\nu }_e$	0.15	0.14	28.57	45.0
$2.54\times {10}^{30}$	$1.22\times {10}^{12}$	${}_{}{}^{48}\mathrm{Si}\to {}_{}{}^{30}O+\Delta N\cdot n-6e^{-}+2{\nu }_e$
		${}_{}{}^{30}O+{}_{}{}^{30}O\to {}_{}{}^{51}\mathrm{Si}+\Delta N\cdot n-2e^{-}+2{\nu }_e$	0.54	0.68	32.64	908.1
$5.78\times {10}^{30}$	$3.73\times {10}^{12}$	${}_{}{}^{53}\mathrm{Si}\to {}_{}{}^{32}O+\Delta N\cdot n-6e^{-}+2{\nu }_e$
		${}_{}{}^{32}O+{}_{}{}^{32}O\to {}_{}{}^{62}S+\Delta N\cdot n$	0.72	0.23	35.47	355.9
$8.69\times {10}^{30}$	$6.16\times {10}^{12}$	${}_{}{}^{64}S\to {}_{}{}^{57}\mathrm{Si}+\Delta N\cdot n-2e^{-}+2{\nu }_e$	0.74	0.03	37.74	3.5
$3.20\times {10}^{31}$	$1.65\times {10}^{13}$	${}_{}{}^{65}\mathrm{Si}\to {}_{}{}^{40}O+\Delta N\cdot n-6e^{-}+2{\nu }_e$
		${}_{}{}^{40}O+{}_{}{}^{40}O\to {}_{}{}^{76}S+\Delta N\cdot n$	0.83	0.05	43.8	98.2
$1.85\times {10}^{32}$	$7.26\times {10}^{13}$	${}_{}{}^{91}S\to {}_{}{}^{86}P+\Delta N\cdot n-e^{-}+{\nu }_e$	0.81	0.006	69.10	0

Table 2. Table 2: Results of our simulations for a parabolic mound profile with a magnetic field B = 10 8 ​ G 𝐵 superscript 10 8 G B=10^{8}\rm G , Accretion rate: 6 × 10 − 13 ​ M ⊙ ​ yr − 1 6 superscript 10 13 subscript 𝑀 direct-product superscript yr 1 6\times 10^{-13}M_{\odot}\;\rm yr^{-1} . The polar spread corresponding to this field strength is 14.9 degrees and the maximum mound height that can be obtained is of 0.2 m, after which the models are unstable. We also compute the ellipticity for this model, although the reader should keep in mind that this is purely the contribution induced by the magnetically supported mound, and does not include contributions due to crustal reactions. The maximum base density is the maximum that is obtained scanning over all grid points of a simulations, while the maximum and minimum of ρ 0 subscript 𝜌 0 \rho_{0} represent the spherical l = 0 𝑙 0 l=0 component of the multipole expansion of the density.

Mound Height	Max base density	Accreted mass	${\rho }_2/{\rho }_0$ (Max)	${\rho }_2/{\rho }_0$ (Min)	${\rho }_0$ (Max)	${\rho }_0$ (Min)	Ellipticity
(m)	($g{\mathrm{cm}}^{-3}$)	($M_{\odot }$)			($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
0.05	$1.808\times {10}^6$	$2.44\times {10}^{-15}$	2.2251	2.1605	$6.323\times {10}^5$	$2.933\times {10}^2$	$4.93\times {10}^{-17}$
0.10	$5.114\times {10}^7$	$1.23\times {10}^{-14}$	2.2249	2.1605	$1.788\times {10}^6$	$8.415\times {10}^2$	$2.44\times {10}^{-16}$
0.15	$9.395\times {10}^7$	$3.48\times {10}^{-14}$	2.2237	2.1605	$3.286\times {10}^6$	$1.685\times {10}^3$	$7.24\times {10}^{-16}$
0.20	$1.446\times {10}^8$	$9.34\times {10}^{-14}$	2.2150	2.1425	$5.104\times {10}^6$	$3.930\times {10}^3$	$2.46\times {10}^{-15}$
0.25	unstable

Table 3. Table 3: Results of the simulation for a magnetic field strength of B = 10 8 𝐵 superscript 10 8 B=10^{8} G and accretion rate M ˙ = 6 × 10 − 13 ​ M ⊙ ​ yr − 1 ˙ 𝑀 6 superscript 10 13 subscript 𝑀 direct-product superscript yr 1 \dot{M}=6\times 10^{-13}M_{\odot}\;\rm yr^{-1} , but a hollow mound profile. All other parameters are set as in table 2 .

Mound Height	Max base density	Accreted mass	${\rho }_2/{\rho }_0$ (Max)	${\rho }_2/{\rho }_0$ (Min)	${\rho }_0$ (Max)	${\rho }_0$ (Min)	Ellipticity
(m)	($g{\mathrm{cm}}^{-3}$)	($M_{\odot }$)			($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
0.040	$1.294\times {10}^7$	$1.16\times {10}^{-15}$	2.1252	2.1249	$4.541\times {10}^5$	$2.122\times {10}^2$	$3.82\times {10}^{-17}$
0.046	$1.596\times {10}^7$	$1.63\times {10}^{-15}$	2.1252	2.1249	$5.601\times {10}^5$	$2.619\times {10}^2$	$5.38\times {10}^{-17}$
0.052	$1.918\times {10}^7$	$2.20\times {10}^{-15}$	2.1252	2.1249	$6.732\times {10}^5$	$3.169\times {10}^2$	$7.26\times {10}^{-17}$
0.058	$2.259\times {10}^7$	$2.88\times {10}^{-15}$	2.1251	2.1248	$7.930\times {10}^5$	$3.770\times {10}^2$	$9.52\times {10}^{-17}$
0.064	$2.618\times {10}^7$	$3.68\times {10}^{-15}$	2.1251	2.1247	$9.191\times {10}^5$	$4.431\times {10}^2$	$1.21\times {10}^{-16}$
0.070	$2.225\times {10}^7$	$4.61\times {10}^{-15}$	2.1251	2.1245	$1.051\times {10}^6$	$5.165\times {10}^2$	$1.52\times {10}^{-16}$
0.076	$3.388\times {10}^7$	$5.67\times {10}^{-15}$	2.1251	2.1242	$1.189\times {10}^6$	$5.991\times {10}^2$	$1.88\times {10}^{-16}$
0.082	$3.797\times {10}^7$	$6.89\times {10}^{-15}$	2.1250	2.1239	$1.333\times {10}^6$	$6.933\times {10}^2$	$2.28\times {10}^{-16}$
0.088	$4.221\times {10}^7$	$8.28\times {10}^{-15}$	2.1249	2.1234	$1.482\times {10}^6$	$8.024\times {10}^2$	$2.75\times {10}^{-16}$
0.094	$4.661\times {10}^7$	$9.87\times {10}^{-15}$	2.1252	2.1228	$1.636\times {10}^6$	1.012	$3.29\times {10}^{-16}$
0.100	$5.114\times {10}^7$	$1.16\times {10}^{-14}$	2.1251	2.1220	$1.795\times {10}^6$	9.016	$3.90\times {10}^{-16}$
0.106	$5.581\times {10}^7$	$1.37\times {10}^{-14}$	2.1251	2.1210	$1.959\times {10}^6$	$2.589\times 10$	$4.61\times {10}^{-16}$
0.112	$6.062\times {10}^7$	$1.61\times {10}^{-14}$	2.1250	2.1196	$2.128\times {10}^6$	$5.332\times 10$	$5.44\times {10}^{-16}$
0.118	$6.555\times {10}^7$	$1.88\times {10}^{-14}$	2.1249	2.1178	$2.301\times {10}^6$	$9.388\times 10$	$6.42\times {10}^{-16}$
0.124	$7.062\times {10}^7$	$2.20\times {10}^{-14}$	2.1249	2.1155	$2.479\times {10}^6$	$1.512\times {10}^2$	$7.58\times {10}^{-16}$
0.130	$7.580\times {10}^7$	$2.57\times {10}^{-14}$	2.1249	2.1124	$2.661\times {10}^6$	$2.303\times {10}^2$	$8.97\times {10}^{-16}$
0.136	unstable

Table 4. Table 4: Results of the simulation for a magnetic field strength of B = 10 8 𝐵 superscript 10 8 B=10^{8} G and accretion rate M ˙ = 3 × 10 − 8 ​ M ⊙ ​ yr − 1 ˙ 𝑀 3 superscript 10 8 subscript 𝑀 direct-product superscript yr 1 \dot{M}=3\times 10^{-8}M_{\odot}\;\rm yr^{-1} , (approximately the Eddington limit), for a parabolic mound profile. The polar spread is now 60.66 degrees at this high accretion rate.

Mound Height	Max base density	Accreted mass	${\rho }_2/{\rho }_0$ (Max)	${\rho }_2/{\rho }_0$ (Min)	${\rho }_0$ (Max)	${\rho }_0$ (Min)	Ellipticity
(m)	($g{\mathrm{cm}}^{-3}$)	($M_{\odot }$)			($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
0.24	$1.901\times {10}^8$	$1.40\times {10}^{-12}$	2.0987	1.2938	$8.9918\times {10}^7$	$3.8860\times {10}^4$	$3.52\times {10}^{-13}$
0.25	$2.022\times {10}^8$	$1.58\times {10}^{-12}$	2.0962	1.2938	$9.5595\times {10}^7$	$4.1953\times {10}^4$	$3.99\times {10}^{-13}$
0.26	$2.144\times {10}^8$	$1.77\times {10}^{-12}$	2.0931	1.2938	$1.0139\times {10}^8$	$4.5307\times {10}^4$	$4.53\times {10}^{-13}$
0.27	$2.269\times {10}^8$	$1.98\times {10}^{-12}$	2.0894	1.2938	$1.0729\times {10}^8$	$4.8980\times {10}^4$	$5.16\times {10}^{-13}$
0.28	$2.396\times {10}^8$	$2.23\times {10}^{-12}$	2.0849	1.2938	$1.1331\times {10}^8$	$5.3043\times {10}^4$	$5.91\times {10}^{-13}$
0.29	$2.526\times {10}^8$	$2.52\times {10}^{-12}$	2.0793	1.2907	$1.1943\times {10}^8$	$5.7600\times {10}^4$	$6.83\times {10}^{-13}$
0.30	$2.657\times {10}^8$	$2.85\times {10}^{-12}$	2.0724	1.2733	$1.2566\times {10}^8$	$6.2793\times {10}^4$	$7.99\times {10}^{-13}$
0.31	$2.791\times {10}^8$	$3.26\times {10}^{-12}$	2.0635	1.2331	$1.3199\times {10}^8$	$6.8832\times {10}^4$	$9.51\times {10}^{-13}$
0.32	unstable

Table 5. Table 5: Results of the simulation for a magnetic field strength of B = 10 10 𝐵 superscript 10 10 B=10^{10} G and accretion rate M ˙ = 3 × 10 − 8 ​ M ⊙ ​ yr − 1 ˙ 𝑀 3 superscript 10 8 subscript 𝑀 direct-product superscript yr 1 \dot{M}=3\times 10^{-8}M_{\odot}\;\rm yr^{-1} , (approximately the Eddington limit), for a parabolic mound profile. The polar spread is 16 degrees.

Mound Height	Max base density	Accreted mass	${\rho }_2/{\rho }_0$ (Max)	${\rho }_2/{\rho }_0$ (Min)	${\rho }_0$ (Max)	${\rho }_0$ (Min)	Ellipticity
(m)	($g{\mathrm{cm}}^{-3}$)	($M_{\odot }$)			($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
1.05	$1.740\times {10}^9$	$4.59\times {10}^{-12}$	2.2315	2.1492	$7.0045\times {10}^7$	$6.0412\times {10}^2$	$1.05\times {10}^{-13}$
1.10	$1.866\times {10}^9$	$5.21\times {10}^{-12}$	2.2308	2.1492	$7.5106\times {10}^7$	$1.0863\times {10}^3$	$1.21\times {10}^{-13}$
1.15	$1.994\times {10}^9$	$5.90\times {10}^{-12}$	2.2302	2.1492	$8.0285\times {10}^7$	$1.7004\times {10}^3$	$1.38\times {10}^{-13}$
1.20	$2.126\times {10}^9$	$6.67\times {10}^{-12}$	2.2297	2.1492	$8.5576\times {10}^7$	$2.4445\times {10}^3$	$1.58\times {10}^{-13}$
1.25	$2.260\times {10}^9$	$7.52\times {10}^{-12}$	2.2292	2.1492	$9.0983\times {10}^7$	$3.3190\times {10}^3$	$1.80\times {10}^{-13}$
1.30	$2.397\times {10}^9$	$8.49\times {10}^{-12}$	2.2287	2.1492	$9.6495\times {10}^7$	$4.3264\times {10}^3$	$2.07\times {10}^{-13}$
1.35	$2.537\times {10}^9$	$9.59\times {10}^{-12}$	2.2283	2.1487	$1.0211\times {10}^8$	$5.4724\times {10}^3$	$2.39\times {10}^{-13}$
1.40	$2.679\times {10}^9$	$1.09\times {10}^{-11}$	2.2341	2.1473	$1.0784\times {10}^8$	$3.6423\times {10}^1$	$2.78\times {10}^{-13}$
1.45	$2.824\times {10}^9$	$1.23\times {10}^{-11}$	2.2330	2.1445	$1.1367\times {10}^8$	$2.0099\times {10}^2$	$3.27\times {10}^{-13}$
1.50	$2.971\times {10}^9$	$1.41\times {10}^{-11}$	2.2321	2.1399	$1.1960\times {10}^8$	$4.9756\times {10}^2$	$3.91\times {10}^{-13}$
1.55	$3.121\times {10}^9$	$1.62\times {10}^{-11}$	2.2314	2.1340	$1.2563\times {10}^8$	$9.2825\times {10}^2$	$4.70\times {10}^{-13}$
1.60	unstable

Table 6. Table 6: Results of the linear extrapolation of the ratio ρ 2 / ρ 0 subscript 𝜌 2 subscript 𝜌 0 \rho_{2}/\rho_{0} at higher densities, following the relation in ( 36 ), ρ 0 = A ​ ρ 2 ρ 0 + B subscript 𝜌 0 𝐴 subscript 𝜌 2 subscript 𝜌 0 𝐵 \rho_{0}=A\frac{\rho_{2}}{\rho_{0}}+B , for the parabolic mound models in table 2 . It can be concluded from the ρ 0 subscript 𝜌 0 \rho_{0} value for ρ 2 / ρ 0 = 1 % subscript 𝜌 2 subscript 𝜌 0 percent 1 \rho_{2}/\rho_{0}=1\% in last column that the ratio ρ 2 / ρ 0 subscript 𝜌 2 subscript 𝜌 0 \rho_{2}/\rho_{0} will be already less than 1 % percent 1 1\% for densities well below 10 9 ​ g ​ cm − 3 superscript 10 9 g superscript cm 3 10^{9}\rm g\;cm^{-3} and thus at densities lower than those of the first deep crustal heating reactions in table 1 .

Accreted mass	Mound Height	A(slope)	B(intercept)	${\rho }_0$ for (${\rho }_2/{\rho }_0=1\%$)
($M_{\odot }$)	(m)	($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
$2.44\times {10}^{-15}$	0.05	$-2.382\times {10}^7$	$5.209\times {10}^7$	$5.2053\times {10}^7$
$1.23\times {10}^{-14}$	0.10	$-7.094\times {10}^7$	$1.550\times {10}^8$	$1.5493\times {10}^8$
$3.48\times {10}^{-14}$	0.15	$-1.857\times {10}^8$	$4.043\times {10}^8$	$4.0402\times {10}^8$
$9.34\times {10}^{-14}$	0.20	Not linear

Table 7. Table 7: Results of the linear extrapolation as in table ( 6 ), but for the hollow mound models in table 3 . The ratio ρ 2 / ρ 0 subscript 𝜌 2 subscript 𝜌 0 \rho_{2}/\rho_{0} remains sizeable at higher densities than for a parabolic mount, up to densities at which reactions in table 1 can occur.

Accreted mass	Mound Height	A(slope)	B(intercept)	${\rho }_0$ for (${\rho }_2/{\rho }_0=1\%$)
($M_{\odot }$)	(m)	($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
$1.16\times {10}^{-15}$	0.040	$-4.111\times {10}^9$	$8.735\times {10}^9$	$8.7321\times {10}^9$
$1.31\times {10}^{-15}$	0.042	$-5.231\times {10}^9$	$1.112\times {10}^{10}$	$1.1113\times {10}^{10}$
$1.46\times {10}^{-15}$	0.044	$-7.136\times {10}^9$	$1.516\times {10}^{10}$	$1.5155\times {10}^{10}$
$1.63\times {10}^{-15}$	0.046	$-9.721\times {10}^9$	$2.065\times {10}^{10}$	$2.0644\times {10}^{10}$
$1.81\times {10}^{-15}$	0.048	Not linear

Table 8. Table 8: Results of the linear extrapolation as in table ( 6 ), for the parabolic mound models in table 4 with accretion rates close to the Eddington limit and B = 10 8 𝐵 superscript 10 8 B=10^{8} G. As can be seen the ratio ρ 2 / ρ 0 subscript 𝜌 2 subscript 𝜌 0 \rho_{2}/\rho_{0} is larger at higher densities than for the weaker accretion rates considered in table ( 6 ), but vanishes just before reaching the densities at which reactions in table 1 can occur.

Accreted mass	Mound Height	A(slope)	B(intercept)	${\rho }_0$ for (${\rho }_2/{\rho }_0=1\%$)
($M_{\odot }$)	(m)	($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)	($g{\mathrm{cm}}^{-3}$)
$1.40\times {10}^{-12}$	0.24	$-3.735\times {10}^8$	$5.715\times {10}^8$	$5.7118\times {10}^8$
$1.58\times {10}^{-12}$	0.25	$-4.381\times {10}^8$	$6.610\times {10}^8$	$6.6047\times {10}^8$
$1.77\times {10}^{-12}$	0.26	$-5.111\times {10}^8$	$3.800\times {10}^8$	$7.5983\times {10}^8$
$1.98\times {10}^{-12}$	0.27	$-6.590\times {10}^8$	$5.176\times {10}^8$	$9.5642\times {10}^8$
$2.23\times {10}^{-12}$	0.28	Not linear

Table 9. Table 9: Results of the linear extrapolation of the ratio ρ 2 / ρ 0 subscript 𝜌 2 subscript 𝜌 0 \rho_{2}/\rho_{0} at higher densities, following the relation in eq.( 36 ), ρ 0 = A ​ ρ 2 ρ 0 + B subscript 𝜌 0 𝐴 subscript 𝜌 2 subscript 𝜌 0 𝐵 \rho_{0}=A\frac{\rho_{2}}{\rho_{0}}+B , for the parabolic mound models in table 5 with accretion rate at the Eddington limit and B = 10 10 𝐵 superscript 10 10 B=10^{10} G.

Accreted mass	Mound Height	A(slope)	B(intercept)	${\rho }_0$ for (${\rho }_2/{\rho }_0=1\%$)
($M_{\odot }$)	(m)			($g{\mathrm{cm}}^{-3}$)
$4.59\times {10}^{-12}$	1.05	$-3.286\times {10}^9$	$7.129\times {10}^9$	$7.1258\times {10}^9$
$5.21\times {10}^{-12}$	1.10	$-3.785\times {10}^9$	$8.207\times {10}^9$	$8.2034\times {10}^9$
$5.90\times {10}^{-12}$	1.15	$-4.437\times {10}^9$	$9.665\times {10}^9$	$9.6605\times {10}^9$
$6.67\times {10}^{-12}$	1.20	$-5.464\times {10}^9$	$1.182\times {10}^{10}$	$1.1821\times {10}^{10}$
$7.52\times {10}^{-12}$	1.25	Not linear

Equations

$$Q_{22}=1.5\times 10^{35}{R}_{10}^4\left(\frac{\delta T_q}{10^5}\right){\left(\frac{E_{th}}{30\mathrm{MeV}}\right)}^3\mathrm{g}\mathrm{c}{\mathrm{m}}^2$$

$$\delta T_q\approx 2.9\times 10^5{R}_{10}^{-4}{I}_{45}{\left(\frac{E_{th}}{30\mathrm{MeV}}\right)}^{-3}\text{mbox}K.$$

$$\Delta T\approx 10^2{C}_k^{-1}{p}_{30}^{-1}{Q}_M\Delta M_{21}\text{mbox}K$$

$$\frac{\delta T_q}{\Delta T}\gtrsim 0.03\text{mbox}(atneutrondrip)$$

$$Q_T\approx 1.5\times 10^{32}\Delta M_{21}{R}_{10}^4\sum _i\frac{(Q_M{)}_i(E_{30}^3{)}_i}{(p_{30}{)}_i(C_k{)}_i}{\left(\frac{\delta T_q}{\Delta T}\right)}_i$$

$${\text{B}}_p=\frac{\nabla \psi (r,\theta )\times \hat{\mathbf{\varphi }}}{r\sin \theta }$$

$$\psi =r\sin \theta A_{\theta }$$

$${\text{j}}_{\varphi }=-\frac{c}{4\pi }\frac{{\Delta }^2\psi }{r\sin \theta }$$

$${\Delta }^2=\frac{{\partial }^2}{\partial r^2}+\frac{\sin \theta }{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\right)$$

$$\nabla p+\rho \nabla \Phi =\frac{\text{j}\times \text{B}}{c}=\rho \nabla F$$

$$\mathbf{g}=-1.86\times 10^{14}\left(\frac{M_s}{1.4M_0}\right){\left(\frac{R_s}{10\mathrm{km}}\right)}^{-2}\mathrm{cm}{\mathrm{s}}^{-2}\hat{\text{r}}$$

$$\frac{1}{\rho }\frac{df}{d\rho }\nabla \rho =\nabla (F(\psi )-\Phi )$$

$$\rho ={\left(\frac{g(\Gamma -1)}{\Gamma K}\right)}^{\frac{1}{\Gamma -1}}{\left(r_0(\psi )-r\right)}^{\frac{1}{\Gamma -1}}$$

$${\Delta }^2\psi =-4\pi r^2{\sin }^2\theta \rho \frac{dF}{d\psi }=-4\pi r^2{\sin }^2\theta \rho g\frac{d{r}_0}{d\psi }$$

$$r_0(\psi )=R_s+r_c\left(1-{\left(\frac{\psi }{{\psi }_A}\right)}^2\right)\text{mbox}parabolicprofile$$

$$r_0(\psi )=R_s+\frac{r_c}{0.25}\left(0.25-{\left(\frac{\psi }{{\psi }_A}-0.5\right)}^2\right)\text{mbox}hollowprofile$$

$$\begin{gathered} r_A=3.53\times 10^3{\left(\frac{B_s}{10^{12}\mathrm{G}}\right)}^{4/7}{\left(\frac{R_s}{10\mathrm{km}}\right)}^{12/7} \\ \times {\left(\frac{\dot{M}}{10^{-9}{M}_{\odot }\mathrm{y}{\mathrm{r}}^{-1}}\right)}^{-2/7}{\left(\frac{M_s}{1.4M_{\odot }}\right)}^{-1/7}\mathrm{km} \end{gathered}$$

$$R_p={\left(\frac{R_s}{r_A}\right)}^{1/2}{R}_s$$

$${\theta }_p={\sin }^{-1}\left(\frac{R_p}{R_s}\right)$$

$$\rho \left(r,cos\theta \right)=\sum _l{\rho }_l(r)Y_{l0}$$

$${\rho }_l(r)=2\pi \int \rho (r,cos\theta )Y_{l0}d\cos \theta$$

$$p(r,\theta )=p_0(r)+\delta p(r,\theta )+O{\left(\delta p\right)}^2$$

$$\delta p(r,\theta )=\frac{\partial p_0}{\partial {\rho }_0}\delta \rho (r,\theta )$$

$$\delta p(r,\theta )=\sum _l\delta p_l(r)Y_{l0}(\theta )$$

$$\delta T_q=-10^2{C}_k^{-1}{p}_{30}^{-1}{Q}_M\Delta M_{21}\left[\left(\frac{\delta p_{30,q}}{p_{30}}\right)+\left(\frac{\delta C_{k,q}}{C_k}\right)\right]$$

$$\frac{\delta T_q}{\Delta T}=-\left[\left(\frac{\delta p_{30,q}}{p_{30}}\right)+\left(\frac{\delta C_{k,q}}{C_k}\right)\right]$$

$$\left(\frac{\delta p_{30,q}}{p_{30}}\right)=\Gamma \frac{\delta {\rho }_q}{\rho }$$

$$C_k=\frac{3}{A}f\left(\frac{T}{{\Theta }_D}\right)$$

$$f(x)=\{\begin{array}{cc}0.8{\pi }^4{x}^3& x\le 0.15\\ 1-0.05x^{-2}& x\ge 0.4\\ 1.7x+0.0083& \text{mbox}otherwise\end{array}$$

$$\delta C_{k,q}=\frac{\partial C_{k,q}}{\partial \rho }\delta {\rho }_q+\frac{\partial C_{k,q}}{\partial T}\delta T_q=\frac{3}{A}\frac{df}{dx}\left(\frac{\delta T_q}{{\Theta }_D}-\frac{x}{2}\frac{\delta {\rho }_q}{\rho }\right)$$

Full text

1 Introduction

The first detection of gravitational waves (Abbott et al., 2016) has opened the field of Gravitational Wave (GW) astronomy and the recent detection of a binary Neutron Star (NS) inspiral (Abbott et al., 2017) has already begun to constrain the Equation of State (EoS) of dense matter (Abbott et al., 2018). NSs are expected to be prolific GW emitters and we not only expect them to play a role in inspirals, but also to emit continuous gravitational waves (CWs), due to unstable modes or non-axisymmetric deformations (dubbed ‘mountains’) that turn the star into a GW pulsar (Riles, 2017).

This last scenario, in particular, has attracted much interest, as on the one side the neutron star crust is thought to be strong enough to support a significant deformation before cracking (Haskell et al., 2006; Horowitz & Kadau, 2009; Baiko & Chugunov, 2018; Johnson-McDaniel & Owen, 2013), on the other accretion onto a magnetised star provides a natural mechanism to produce crustal asymmetries and source a deformation (Bildsten, 1998; Ushomirsky et al., 2000). It has even been suggested that the spin rate of accreting NSs in Low Mass X-ray Binaries (LMXBs) may be set by the competition between the accretion torque that is spinning-up the star and the GW torque that removes angular momentum from the system (Papaloizou & Pringle, 1978). This would explain why we do not observe sub-millisecond pulsars, despite the Keplerian breakup frequency of a NS $\nu_{k}$, being above $\nu_{k}\approx 1200$ Hz for any physically realistic model that is causal in the core and includes a crust (Haskell et al., 2018).

A recent analysis of the spin-distribution of accreting NSs in LMXBs by Patruno et al. (2017) has also revealed that the distribution is bimodal, with a broad population appearing at lower frequencies (and consistent with the observed population of millisecond radio pulsars, of which the LMXBs are thought to be the progenitors Alpar et al. (1982)) and a much more narrowly distributed population of rapidly rotating neutron stars with $550$ Hz $\lesssim\nu\lesssim 720$ Hz. This feature, in particular, suggests a torque that scales strongly with frequency and is difficult to explain simply involving accretion torques, but can be naturally explained if GW torques are acting on the system (Gittins & Andersson, 2018).

One of the pulsars in the fast population, PSR J1023+0038, has attracted particular attention. This is a transitional pulsar, that transitions from being accretion powered in outburst to radio emission during quiescence (Archibald et al., 2009). It has been possible to time this pulsar both in the X-ray and radio phases, and measurements of the spin-down rate have revealed that the NS is spinning down $\approx 27\%$ faster when it’s accreting (Jaodand et al., 2016). This result is somewhat unexpected, as standard accretion models can predict spin-down during accretion, but generally require some fine-tuning for the spin-down rate in both phases to differ so little (Jaodand et al., 2016). Gravitational wave emission, however, offers a scenario that can explain the slight increase in spin-down rate quite naturally if a mountain is being built on the crust during the accretion outburst (Haskell & Patruno, 2017).

If light elements are accreted at the surface, they are subsequently pushed deeper into the crust when new material is added, and here they undergo several reactions, including pycno-nuclear reactions and electron capture, that change the composition and release heat locally (Haensel & Zdunik, 1990). This deep-crustal heating is responsible for heating X-ray transients, which then cool during quiescence as has been observed in many systems (Wijnands et al., 2017). If part of the reaction layers are not exactly axisymmetric, temperature and composition gradients can source a ‘mountain’ and lead to GW emission (Bildsten, 1998; Ushomirsky et al., 2000).

The main parameters that determines the size of the mountain are the total accreted mass, which depend on the accretion rate and outburst duration, and the amount of asymmetry in the heat release. For example Haskell & Patruno (2017) found that if $\approx 2\%$ of the heat released is quadrupolar, then a large enough mountain, that explains the observed spindown rate, can be built on PSR J1023+0038.

How much of the heat emission will be quadrupolar, or how quadrupolar the surfaces of equal composition will be, is however, generally unknown. One can put an upper limit of around $50\%$ on the ratio between quadrupolar and spherical components of the heat deposition from the fact that no pulsations are visible in quiescence with an amplitude larger than $\approx 50\%$, while if there were quadrupolar heat deposition, the heat would diffuse out in quiescence and cause pulsations at twice the rotation frequency with similar amplitude (see Haskell et al. (2015) for a detailed analysis). There is, ultimately, no realistic estimate of how large the quardupolar component is expected to be, thus limiting the predictive power of this model.

In this paper we consider a natural source of asymmetry, namely the magnetic field. We model accretion onto the polar cap of a neutron star numerically with the code of Mukherjee & Bhattacharya (2012), and calculate how distorted the density and pressure profiles will be in the outer layers of the star. We then extrapolate to higher densities, to calculate the quadrupolar component of the temperature perturbation due to reactions in the outer crust up to neutron drip density ($\rho\approx 10^{11}$ g cm*-3*). From this we use the results of Ushomirsky et al. (2000) to estimate the induced ellipticity and gravitational wave signal. In general we find that only relatively small ellipticities can be produced by this mechanism, smaller than could be detected by Advanced LIGO and Virgo and smaller than what would be required to explain the additional spindown of PSR J1023+0038.

2 Deep crustal heating

A NS in an LMXB accretes matter, generally hydrogen or helium, from the companion, and as accretion proceeds these light nuclei are pushed down to higher densities, and can undergo a number of reactions, including electron captures and pycno-nuclear fusions, which produce heavier nuclei and push the composition towards equilibrium (Haensel & Zdunik, 1990). This process releases a few MeV per accreted baryon of heat in the crust, which is thought to power the quiescent luminosity of NS in LMXBs and can explain the observed cooling of transient LMXBs (Wijnands et al., 2017). If the reaction layers are not perfectly axisymmetric, gradients in temperature and composition resulting from asymmetric accretion can source a quadrupolar deformation, sustained by the rigidity of the crust (Bildsten, 1998). Ushomirsky et al. (2000) examined this problem in detail and calculated the quadrupole induced by reactions in the outer crust. Their results, assuming a quadrupolar temperature increase $\delta T_{q}$, can be approximated as:

[TABLE]

where $R_{10}$ is the radius in units of 10 km, $E_{th}=\mu_{e}$, with $\mu_{e}$ the electron chemical potential, is the threshold energy for the pycnonuclear reactions responsible for deep crustal heating in a given layer. Since the required quadrupole $Q_{22}$ to explain the behaviour of PSR J1023+0038 is $4.4\times 10^{35}I_{45}$ g $\rm{cm^{2}}$ Haskell et al. (2018), where $I_{45}$ is the moment of inertia of the NS in units of $10^{45}$ g $\rm{cm^{2}}$, from eq. (1) we see that the required quadrupolar temperature variation for this star is

[TABLE]

Assuming that deep crustal heating is the only source of temperature increase, the total increase in temperature is given as Ushomirsky & Rutledge (2001):

[TABLE]

where $C_{k}$ is the heat capacity in units of the Boltzman constant per baryon, $p_{30}$ is the pressure in units of $10^{30}$ erg/$\rm{cm^{3}}$, $Q_{M}$ is the heat released locally by the reactions per accreted baryon, in units of MeV, and $M_{21}$ is the accreted mass in units of $10^{21}$ g.

To obtain an estimate for PSR J1023+0038 Haskell & Patruno (2017) considered a simplified model with only the contribution due to the most energetic reactions close to the neutron drip layer, at $\rho\approx 4\times 10^{11}$ g cm*-3*, although note that for accreted crusts this may be shifted to higher densities depending on the composition of the ashes of X-ray bursts, as discussed by Chamel et al. 2015. We have used $\Delta M_{21}=0.1$ following Haskell & Patruno (2017), which is reasonable for a month of accretion. In this case, taking $E_{th}=30$ MeV, $Q_{M}=0.5$ MeV, $C_{k}=10^{-6}$ and $p_{30}=1$, one obtains

[TABLE]

Naturally a more detailed analysis should consider the contribution due to all reaction layers, which will add up to contribute to the total quadrupole, thus reducing the estimate in (4). Taking the results of Haskell & Patruno (2017) and summing the contribution of all layers in the outer crust, as suggested by Ushomirsky et al. (2000), for a linearized model, the total quadrupole generated by $i$ reactions in different layers will be given as:

[TABLE]

where $E_{30}$ is $E_{th}$ in units of $30$ MeV, and all quantities with the subscript $i$ must be calculated at the pressure and density corresponding to the reaction layer $i$.

It remains an open issue whether deformations of order $\delta T_{q}/\Delta T\approx 0.01$ can occur in the crust of an accreting NS, and without an estimate of this quantity it is impossible to obtain a meaningful theoretical estimate of the quadrupole from (5), let alone compare this to observations.

There is, however, a natural source of asymmetry in the system, as the NS is magnetised, and unless the field is too weak, matter is accreted onto the polar caps and then spreads due to lateral pressure gradients, that are opposed by magnetic stresses. This can lead to the creation of magnetic mounds in accreting systems (Payne & Melatos, 2004, 2007; Mukherjee, 2017), and significantly deform the field structure, leading to burial of the dipolar component (Shibazaki et al., 1989). In fact, for high magnetic fields and accretion rates, the magnetically confined mountain may itself be large enough to lead to a detectable gravitational signal (Priymak et al., 2011; Haskell et al., 2015). We will see, however, that for the weaker fields we consider the magnetic contribution to the ellipticity is generally much smaller than the ‘thermal’ contribution due to deep crustal heating.

In summary our strategy consists of splitting the full problem, which would include contributions due to magnetic stresses, elasticity and the readjustment of the reaction layers due to the non-barotropic nature of a realistic EoS in the crust, into following steps. First we calculate barotropic magnetic equilibrium to estimate the temperature asymmetries in the crust of the accreting star. Then we use these to calculate the source terms that lead to the readjustment of the capture layers and then using the results of Ushomirsky et al. (2000), we obtain the elastic response of the crust.

3 Temperature asymmetries

To calculate the temperature asymmetries in the crust we consider a sequence of static equilibria of a mound of accreted matter confined by a magnetic field in the neutron star outer layers. The equilibrium solution is constructed by numerically solving the Grad-Shafranov (GS) equation using the numerical GS solver developed by Mukherjee & Bhattacharya (2012), described in detail in Section 3.1. As in the original algorithm of the GS solver, the effects of continued accretion are not taken into account. Although it has been demonstrated that pressure driven MHD instabilities Mukherjee et al. (2013a, b) can occur beyond a threshold mass of the magnetically confined mounds, Vigelius & Melatos (2008) have shown that the instabilities may saturate to a new equilibrium configuration. For analytical simplicity, we perform the analysis in this paper on the equilibrium solutions obtained from the GS solver, without exploring the stability of the mound, to obtain an upper limit on the deformations and the resultant gravitational wave signals. We then, use the equilibrium solution to estimate the asymmetries in the crust which is described in Section 3.2

3.1 Formulation of Grad-Shafranov equation

To simulate the mass and magnetic field configuration on the surface of a neutron star in an LMXB we consider the equilibrium profile of a ‘mountain’ on the surface of a star accreting matter from a disc, where the magnetosphere cuts off the accretion disk at the Alfven radius, the distance from the star where the energy density of the stellar magnetic field is balanced with the energy density of the inflowing matter and the matter is accreted onto the poles. As we shall see in the following this allows us to define the extent of the polar cap on which matter is deposited by accretion. More detailed treatments of the inner disc radius are possible (Ghosh & Lamb, 1979; Wang, 1996; Spruit & Taam, 1993; Andersson et al., 2005; Kluźniak & Rappaport, 2007), but generally the position of the inner edge of the disc $r_{inner}$ is found to be proportional to the Alfven radius $r_{A}$, such that $r_{inner}=\xi r_{A}$, with $0.4\lesssim\xi\lesssim 1$ (D’Angelo, 2017). For simplicity, and given the other sources of uncertainty in our treatment, we will thus simply take $r_{inner}=r_{A}$.

The initial magnetic field is dipolar and the polar caps are assumed to be axisymmetric about the magnetic axis z. For zero toroidal magnetic field in an axisymmetric system $\left(r,\theta,\phi\right)$, in ideal magnetohydrodynamics (MHD) the poloidal magnetic field is given as:

[TABLE]

where $\psi$, the flux function is given as:

[TABLE]

where, $A_{\theta}$ is the azimuthal component of the vector potential and the system is symmetric with respect to azimuthal coordinate $\phi$. The azimuthal component of the current is then,

[TABLE]

[TABLE]

where, $\Delta^{2}$ is the Grad-Shafranov operator in spherical coordinates. For a system in static equilibrium, the Euler equation is,

[TABLE]

where $\Phi$ is the gravitational potential and $\boldsymbol{\nabla}F=\frac{\boldsymbol{\nabla}p}{\rho}+\boldsymbol{\nabla}\Phi$, so F is a function constant along a flux surface. For a spherical coordinate system and Newtonian gravity with constant acceleration,

[TABLE]

assuming $p=f(\rho)$, the density profile is obtained by integrating eq.(10) to get,

[TABLE]

For a polytropic equation of state $p=K\rho^{\Gamma}$ eq.(12) gives,

[TABLE]

where $r_{0}=F(\psi)/g$ is an arbitrary mound height function, the choice of which gives the shape of the mound. The constant of integration is taken to be zero. Substituting eq.(13) in eq.(10) gives,

[TABLE]

Eq.(14) is the Grad-Shafranov equation in spherical coordinates. The above equation has been solved using a parabolic mound height function of the form:

[TABLE]

, and hollow mound height function given as:

[TABLE]

which accounts for the fact that mass loading will occur over a finite range of accretion disc radii, as discussed by Mukherjee et al. (2013a). Here $\psi_{A}$ is the flux function at the Alfven radius and $r_{c}$ the cutoff chosen for the mound height. The density is determined by mound height profile and eq.(13).

Specifically we take polytropic equation of state (EOS) of the form : $p(\rho)=K\rho^{\Gamma}$ where $\Gamma$ is the adiabatic index and $K$ is measured in cgs units (dyn $\rm g^{-\Gamma}\rm cm^{3\Gamma-2}$). The values of $K=5.4\times 10^{9}$ and $\Gamma=5/3$ are chosen to crudely approximate to the density regimes of interest in the crust i.e degenerate neutron gas $10^{12}\leq\rho/$(g $\rm cm^{-3})\leq 10^{16}$. The neutron star parameters are taken as $M_{s}=1.4M_{\odot}$, $R_{s}=10\rm km$, where $M_{s}$ and $R_{s}$ are the mass and radius of the neutron star and $M_{\odot}$ is the solar mass. We note that our choice of EOS is mainly dictated by the desire to simulate high densities and models with large values of the accreted mass. Our choice ensures this (see Priymak et al. (2011) for a more detailed discussion) but is not entirely consistent as, in practice, models used for our extrapolation do not reach such high densities. A consistent model should consider the pressure as mainly due to degenerate electrons rather than neutrons. We computed equilibria using various EOS and found that a degenerate electron EOS is suitable for modelling the accreted matter only in the upper layers and for small accreted mass. To explore the maximum limit of the accreted mass sustainable inside the crust we use the non-relativistic degenerate neutron EOS to model the possible confinement of denser matter, keeping in mind that it is likely to provide an upper limit to the true value of the quadrupole that we calculate, and that a more realistic EOS should be considered in future work.

The Alfven radius $r_{A}$ is deduced by equating ram pressure to magnetic pressure and can be expressed as Elsner & Lamb (1977)

[TABLE]

The maximum size of the polar cap is determined by the field lines from the Alfven radius on to the neutron star surface. The polar cap radius $R_{p}$ can be expressed in terms of Alfven radius $r_{A}$ as :

[TABLE]

If $R_{p},\;\theta_{p}$ are the polar cap radius and the opening polar cap angle respectively, then $\theta_{p}$ is given as:

[TABLE]

3.2 Estimate of asymmetries in the crust

In the spherical coordinate system $\left(r,\theta,\phi\right)$ used here, with $r$ the distance from the center of the star, $\theta$ the angle from the magnetic axis, and $\phi$ the azimuthal angle, the system is symmetric with respect to the azimuthal coordinate $\phi$. The density values can be written as $\rho\left(r,cos\theta\right)$, and expanded in terms of standard spherical harmonics with $m=0$ ($Y_{l0}$) in the form:

[TABLE]

and each component is calculated as :

[TABLE]

As we are interested in the quadrupolar deformation which gives the leading order contribution to the gravitational wave flux, we focus on the $l=2$ coefficients and calculate the values for the ratio, $\rho_{2}/\rho_{0}$ for each layer in the mound. This quadrupolar variation in density generated due to magnetic stress is then used to evaluate quadrupolar variation in the temperature for the estimation of the quantity $\delta T_{q}/\Delta T$, from equations (26), (27) and (30).

We note that as the numerical setup of the system is symmetric with respect to the azimuthal coordinate $\phi$, a direct expansion in Legendre polynomials may have been more natural. However we choose to expand in spherical harmonics as the true quantity that is needed to estimate the quadrupolar deformation in eq. (2) is the $l=2,m=2$ component of the temperature perturbation. We make the assumption as in Haskell et al. (2008) that $\delta T_{22}\approx\delta T_{20},\,\delta\rho_{22}\approx\delta\rho_{20}$, the subscripts being the $l$ and $m$ values respectively. We have chosen to do so since for slow rotation of the star the two quantities differ only by a geometric factor due to the misalignment of the rotation axis (which sets the spherical coordinate basis with respect to which the spherical harmonics are defined) and the magnetic axis. The estimate made here thus gives the $l=2,m=0$ component of the temperature perturbation rather than the $l=2,m=2$ component, but allows us to solve for more tractable model in which the magnetic and rotational axis are aligned (Bonazzola & Gourgoulhon, 1996; Haskell et al., 2008).

Keeping in mind the approximation $\delta T_{22}\approx\delta T_{20}$, to estimate the quadrupolar perturbation in the temperature, let us consider linear perturbations of the pressure, around a spherically symmetric background value $p_{0}(r)$ i.e,

[TABLE]

where $\delta p(r,\theta)$ can also be written in terms of density perturbations as,

[TABLE]

If we assume an expansion in spherical harmonics for the perturbations also, such that , e.g.

[TABLE]

we can use this linear perturbation in eq. (3) to obtain the value for quadrupolar temperature variation as,

[TABLE]

or, equivalently:

[TABLE]

where with the subscript $q$ we are now indicating the $l=2$ coefficient of the expansion. We assume $Q_{M}$ to be constant which is a good approximation on the timescale of an accretion outburst as $Q_{M}$ can only change on much longer evolutionary timescales when the star first starts accreting and the cold catalyzed crust is replaced by accreted material. Assuming a polytropic equation of state we can write

[TABLE]

At the densities of interest the main contribution to the heat capacity is that of the ions in the lattice (Potekhin et al., 2015). We therefore take (Chong & Cheng, 1994):

[TABLE]

where $\Theta_{D}=3.48\times 10^{3}\rho^{1/2}Z/A$ K is the Debye temperature at a given density $\rho$, is function of the proton number $Z$ and the atomic number $A$. The Debye function $f(x)$, with $x=T/\Theta_{D}$, can be approximated as (van Riper, 1991)

[TABLE]

Perturbing eq. (28) we have,

[TABLE]

so that for the temperatures of interest i.e. $T\lesssim 10^{7}$K or $x\leq 0.4$ we have

[TABLE]

In the following we will assume that the crustal heating is strong enough in the capture layers, so that $T\approx\Delta T(l=0)$ (Ushomirsky & Rutledge, 2001). We obtain $\frac{\delta T_{q}}{\Delta T}$ by combining eq. 31, 26 and 27. Thus the spherically symmetric heat increase obtained from equation (3) sets the background temperature of our model, and also identify the background $\rho$ with $\rho(l=0)$ extracted from our simulation. Thus the main assumption that allows us to obtain the result is that ratio $\delta T_{q}/\Delta T$ can be expressed in terms of the ratio of $l=2$ and $l=0$ component of mass density perturbation $\rho_{2}/\rho_{0}$, which is generated due to the confinement of accreted material by magnetic stress produced by the bending of field lines. It is based on the fact that the temperature of the crust in many cases is set by crustal heating (Ushomirsky & Rutledge, 2001). The only additional ingredient we need for our calculation of the quadrupole is the density $\rho$, pressure $p$, threshold energy $E_{th}$ and energy release $Q_{M}$ of the deep crustal heating reactions in the outer crust. We take the recent results of Fantina et al. (2018), and the relevant parameters can be found in table 1.

4 Results

To determine the quadrupolar density deformations we computed equilibria for varying values of accreted mass $\Delta M$ and magnetic field strength $B$, for two different values of the accretion rate, a low value of $\dot{M}=6\times 10^{-13}M_{\odot}\;\rm yr^{-1}$ which would be appropriate for a system such as J1023+0038 in which continuous gravitational wave emission was suggested (Haskell & Patruno, 2017), and a high value of $\dot{M}=3\times 10^{-8}M_{\odot}\;\rm yr^{-1}$, more appropriate for persistent sources accreting close to the Eddington limit, which are likely to be the best targets for gravitational wave searches (Haskell et al., 2015).

We compute the density profile for each model and expand it as equation 20, after which we compute the ratio $\rho_{2}/\rho_{0}$. Assuming that the effect of the magnetic field is weak compared to gravity, so as to maintain a spherically symmetric background density, we make the approximation that $\delta\rho_{q}/\rho_{0}\approx\rho_{2}/\rho_{0}$, which allows us to determine $\delta T_{q}/\Delta T$ from equation (26).

In tables 2, 3 and 4 we show example of the results for an initial magnetic field strength of $B=10^{8}$ G (appropriate for J1023+0038) for both the parabolic (filled) and hollow profile. Examples of the ratio $\rho_{2}/\rho_{0}$ are plotted in figures 1 and 2, for different values of the accreted mass (or alternatively mound height), and extensive tests have been run for different values of the magnetic field strength. Table 5 shows the results of the simulation for a magnetic field strength of $B=10^{10}$ G and accretion rate $\dot{M}=3\times 10^{-8}M_{\odot}\;\rm yr^{-1}$, (approximately the Eddington limit), for a parabolic mound profile.

We see that, in general, the densities at the base of the crust are lower than that of many of the reaction layers we consider, and the ratio $\rho_{2}/\rho_{0}$ appears to be linearly decreasing with density. This suggests that we may extrapolate our results to higher densities by means of a fit of the dependence of $\rho_{2}/\rho$ on $\rho$. We choose a linear fit, mostly for simplicity, as this appears to be a good approximation to our data. In fact we do not have enough data points close to the base of the mountain and spread over a large enough dynamical range to meaningfully perform a logarithmic fit, or fit a power-law with varying index. We can further assess this approximation with a simple analysis of the forces acting on the system. Let us consider a plane parallel section of a neutron star, where gravity $g$ is directed in the negative $z$ direction. The Alfven wavelength, i.e. the length-scale at which magnetic stresses become the dominant force in the system Heng & Spitkovsky (2009)

[TABLE]

with $B_{z}$ the magnetic field along the $z$ axis and the scale-height $H=p/\rho g$. We can now estimate the quadrupolar mass distribution by assuming that the accreted material can only spread within $\lambda_{a}$ and make the approximation that $\lambda_{a}$ can be taken to be also the curvilinear length-scale over which matter is confined at a given radius (Cavecchi et al., 2011). Then the $l=2$ component of the density is given by

[TABLE]

which as expected vanishes for $\lambda_{a}/R=0$ and $\lambda_{a}/R=\pi/2$, i.e. when matter has spread to the whole surface. If we expand close to $\lambda_{a}/R=\pi/2$ we see that for high densities $\rho_{2}/\rho\propto(\pi/2-\lambda_{a}/R)$. Assuming a polytropic equation of state, we have from (32)

[TABLE]

which gives $\frac{\rho_{2}}{\rho}\propto\rho^{3/2}$ for our $\Gamma=5/3$ equation of state, thus indicating a powerlaw behaviour that is not, however, too far from our linear approximation. Note that the above analysis is valid for weak fields in which matter can spread significantly. An approximate estimate of the confinement radius can be found by equating the magnetic Lorentz force to the pressure gradient, as derived in Brown & Bildsten (1998)

[TABLE]

Repeating the same analysis as above, we find that for strong fields we may expect $\rho_{2}/\rho\propto\rho^{(2\Gamma-1)}$, which for our polytrope gives $\rho_{2}/\rho\propto\rho^{7/3}$. In practice in the presence of significant bending of the field lines due to large accreted masses, the ratio $\rho_{2}/\rho_{0}$ obtained from our simulations deviates significantly from a linear relation, as can be seen in figures 1 and 2, and the effects of the lower boundary play a significant role. In these cases, we do not attempt to extrapolate our results to higher densities, but rather extrapolate only when the linear fit is justified by the data. An example of a linear fit for parabolic mound profile is shown in figure 3

Overall, our main justification for a linear fit is, thus, numerical. Although a power-law fit of the form $\rho_{2}/\rho=\rho_{0}^{k}$ may be more justified, the analysis above suggests that $1\lesssim k\lesssim 2$, and is generally valid only for small mounds in which the field lines are not bent significantly and matter does spread far from the polar caps. Given the large physical uncertainties in the model for larger mounds, and small dynamical range of the data we can fit over, we choose a linear fit as suggested by the data. If, indeed, the deformation $\rho_{2}/\rho$ scales with a somewhat higher power of density, closer to $k\approx 2$, then our results are to be interpreted as upper limits.

We also test the validity of our extrapolation procedure and of the perturbative expansion by plotting also the ratio between higher $l$ multipoles of density and the spherical component, i.e. the ratios $\rho_{l}/\rho_{0}$, up to $l=10$. The results are given in figure 5 for both the parabolic and hollow mound profiles. We see that for a parabolic mound, in general higher multipoles decrease faster as density increases, validating our approximation that at higher densities at which the reactions in table 1 occur, the density distribution is mainly spherical with small deviations from axisymmetry. For a hollow mound profile the situation is different, as the density distribution is still highly asymmetric at the base of our grid. This is to be expected given the ’hollow’ nature of the profile that naturally leads to a roughly bimodal density distribution. This can be seen clearly by comparing the density distributions in figure 6, where it is clear that even for small accreted masses and mound heights of 0.1 m for a parabolic mound, the density distribution has a single peak centered around the axis, while for a hollow mound the peak is off center, and the density distribution is highly non-spherical. As we will see below this is not only a more realistic setup, but also the most interesting for gravitational wave emission as it naturally leads to higher asymmetries and ‘mountains’. It is also, however, the setup for which our extrapolation is less reliable and our results are thus likely to only be lower limits, as they neglect non-linear effects and the contribution from higher values of $l$. A consistent, non-linear analysis of the problem is clearly needed in the future to assess this problem quantitatively.

The results for the extrapolation are given in tables 6, 7 and 8 where we show the results of our fitting procedure, for which we assume a linear dependence between $\rho_{0}$ and the ratio $\rho_{2}/\rho_{0}$ of the form:

[TABLE]

In general we see that for small amounts of accreted mass, comparable to what a weakly magnetised and slowly accreting source such as PSR J1023+0038 would accrete in a few days of outburst, the quadrupolar deformation, for filled accretion mounds, vanishes at densities of a few times $\rho\approx 10^{8}\rm g\;cm^{-3}$. Deformations are present at higher densities for higher magnetic fields, which limit the spread of matter on the star, but generally are not present above $\rho\approx 10^{9}\rm g\;cm^{-3}$. The same is true also for sources accreting close to the Eddington limit, which suggests that even in persistently accreting sources, which are prime targets for gravitational wave detection (Haskell et al., 2015), one cannot obtain large thermal mountains.

For the more realistic case of a hollow mound, however, deformations can persist up to $\rho\approx 10^{10}\rm g\;cm^{-3}$ in weakly magnetised neutron stars. For strong magnetic fields, however, the the ratio $\rho_{2}/\rho_{0}$ is a strongly non-linear function of $\rho$, as is also the case for large accreted masses, for which the mound height approaches it’s critical value. In all cases where the extrapolation to high densities is possible, i.e. for weak fields and low accretion rates (as is the case for PSR J1023+0038) deformations of the order of $\delta T_{q}/\Delta T\approx 0.01$ are not possible at densities $\rho\approx 10^{12}\rm g\;cm^{-3}$, such as those of the most energetic pycno-nuclear reactions considered by Haskell & Patruno (2017).

5 Induced quadrupole

Following the results of the previous section, it is clear that in most cases no significant deformation persists up to densities at which reactions can occur. For parabolic mounds this is only the case for high accretion rates and strong magnetic fields. For weaker fields, of the order of $B\approx 10^{8}$ G, as commonly inferred in LMXBs, only the hollow mounds allow significant deformation at the densities at which the reactions in table 1 take place.

We consider this case first, and calculate the quadrupole for such a hollow mound. Despite not releasing as much energy as reactions in deeper layers, it is clear that shallow layers play a role, as deformations can be large at those depths. As an example let us consider the case of PSR J1023+0038, for which the increase in spin-down rate during outburst can be explained by a quadrupole of Haskell & Patruno (2017)

[TABLE]

where we assume a moment of inertia $I=10^{45}$ g cm2. We can compare this to the quadrupole of our model with an elliptical mound, $B=10^{8}$ G and an accreted mass of $4\times 10^{18}$ g (approximately the amount of mass PSR J1023+0038 would accrete in the first day or two of outburst).

The maximum density at the base of our grid is $\rho_{0}=5.601\times 10^{5}\rm g\;cm^{-3}$ and the corresponding $\rho_{2}/\rho_{0}=2.1249$. The relation between $\rho_{2}/\rho_{0}$ and $\rho$ appears linear, as can be seen in figure 4 so we extrapolate our results with the linear fit $\rho_{0}=-9.721\times 10^{9}(\rho_{2}/\rho_{0})+2.065\times 10^{10}\rm g\;cm^{-3}$. From this, and using the data from table 1 in equation (5), we obtain, assuming a background temperature $T=10^{7}$ K, $A=85$ and $Z/A=0.4$.

[TABLE]

In comparison the mass quadrupole induced by the matter confined by the magnetic field on the surface is $Q\approx 10^{29}$ g cm2 and negligible. We see that the quadrupole is well below what is needed to explain the observations. One may object that a larger amount of mass is accreted during an outburst, $\Delta M\approx 10^{20}$ g, and despite our models showing that the mound is unstable for such values, numerical time-evolutions suggest that the mound height saturates even if additional mass is added (Vigelius & Melatos, 2008). We may thus extrapolate to larger values of the accreted mass. $Q_{T}$ grows linearly with accreted mass, so even for accreted masses of order $10^{21}$ g, which is generally more than accreted during an outburst the quadrupole would not be large enough to explain the spin-down of PSR J1023+0038.

Let us now turn our attention to the high accretion rate, where the star is accreting close to the Eddington limit and the background field is $B=10^{10}$ G. In this case density perturbations extend to the first reaction layer in table 1. If we assume in this case, a hotter star with an internal temperature of $T=10^{8}$ K, from the results in table 9 we obtain a quadrupole of

[TABLE]

which corresponds to an ellipticity of $\epsilon\approx 5\times 10^{-11}$ and a signal that is likely to be too weak for detection with current ground based interferometers (Abbott et al., 2019a, b, c).

We note, however, that several studies of cooling X-ray transients have found that the cooling curves cannot be reproduced unless a shallow heating source with $Q_{M}\approx 1-10$ MeV is included at densities $\rho\approx 10^{8}-10^{10}\rm g\;cm^{-3}$ (Deibel et al., 2015; Waterhouse et al., 2016; Parikh et al., 2017, 2018). If an additional shallow heating layer is added to the reactions in table 1, at $\rho=10^{9}\rm g\;cm^{-3}$ and with $Q_{M}\approx 5$ MeV, using the BSk20 EoS (Goriely et al., 2010), and take an accreted mass of $\Delta M\approx 10^{20}$ g, we obtain $Q_{T}\approx-8\times 10^{34}$ g cm2. A strong shallow heating source is thus consistent with a scenario in which gravitational waves spin down PSR J1023+0038 during outburst. Furthermore if such hypothetical shallow heating sources are present in systems accreting close to the Eddington limit, with a magnetic field of $B=10^{10}$ G, then one would have $Q_{T}\approx-3\times 10^{38}$ g cm2, corresponding to an ellipticity of $\epsilon\approx 4\times 10^{-7}$ and potentially detectable in the near future by Advanced LIGO and Virgo from known persistently accreting neutron stars in LMXBs.

The above is an estimate of the quadrupole that can be formed during a single outburst. If the quadrupole is indeed sourced uniquely by temperature asymmetries, we expect it to be washed out on a thermal timescale for the crust, which for a capture layer at a pressure $P_{30}$ in units of $10^{30}$ erg/cm3, is given by Brown et al. (1998)

[TABLE]

so that a shallow ‘thermal’ mountain as we have considered would be washed out entirely a few months after the outburst. However, as a consequence of the outburst, compositional asymmetries will build up as the reactions proceed, and may be accumulated over several outbursts (Ushomirsky et al., 2000), allowing to build up a mountain over several outbursts. In fact, the recent analysis of the spin-distribution of accreting neutron star by Gittins & Andersson (2018) shows that population synthesis models including a quadrupole from a mountain built up over many outbursts fit the data better than the models which assume that a mountain is built up over a single outburst and then washed away.

6 Conclusions

We have calculated numerical models of the outer layers of an accreting magnetised neutron star, using the numerical Grad-Shafranov code of Mukherjee & Bhattacharya (2012). The density profile in these models has been expanded in spherical harmonics and the ratio between the quadrupolar and spherical density perturbations $\delta\rho_{2}/\delta\rho_{0}$, extrapolated linearly to higher densities. We find that for a filled mound with a parabolic profile function for the height, the quadrupolar density perturbation vanishes at densities lower than $\rho\approx 10^{9}\rm g\;cm^{-3}$ for field strengths of $B\approx 10^{8}$ G, typical of neutron stars in LMXBs, for all accretion rates, and only for higher magnetic fields of $B=10^{10}$ G do we have significant perturbations up to densities of $\rho\approx 10^{10}\rm g\;cm^{-3}$ . For a more realistic hollow accretion profile (Mukherjee et al., 2013a) the quadrupolar density perturbation persists up to densities at which deep crustal heating reactions occur (Fantina et al., 2018). For these models we calculate the quadrupolar temperature perturbations due to the reactions, following the formulation of Ushomirsky et al. (2000), to estimate the size of the additional mass quadrupole that would be induced. We consider in particular the case of PSR J1023+0038, for which it has been suggested that a gravitational wave mountain may be created during an accretion outburst (Haskell & Patruno, 2017), explaining the enhanced spin-down of the star. For standard reactions we obtain quadrupoles of the order

[TABLE]

for this source, which is not sufficient to explain the additional spin-down.

However the cooling curves of many accreting sources reveal the need to add additional shallow heating sources at densities $\rho\lesssim 10^{10}\rm g\;cm^{-3}$ (Deibel et al., 2015; Waterhouse et al., 2016; Parikh et al., 2017, 2018). Indeed if we add a strong shallow heating source which releases $Q=5$ MeV of heat at a density of $\rho=10^{9}\rm g\;cm^{-3}$, we find that one can build, over the outburst, a mountain corresponding to a quadrupole

[TABLE]

which would explain the increased spin-down of PSR J1023+0038 during outburst. Furthermore, for a system accreting persistently at the Eddington rate and with a magnetic field of $B=10^{10}$ G, the quadrupole could be as large as

[TABLE]

and lead to a GW signal that is potentially detectable by Advanced LIGO and Virgo for known persistently accreting LMXBs, although even at this level the search would be challenging (Watts et al., 2008).

We stress, however, that our estimates rely on extrapolations at densities higher than those of our computational domain, and that the large quadrupoles we obtain, larger than those induced by the magnetic deformations themselves, indicate that treating this effect as a linear perturbation is likely to be an inadequate approach. By limiting ourselves to cases in which the relation between the quadrupolar and spherical density perturbations is linear, we are limiting ourselves to cases in which there is no substantial line bending, and the magnetic deformations are also small. For larger deformations, MHD instabilities may be present (Litwin et al., 2001; Mukherjee et al., 2013a, b) and may limit the growth of mounds beyond the heights considered in this work. In reality the true values of the mass quadrupole may thus be close to the lower limits we compute. Furthermore, the equation of state we use does not consistently model all the regions of the crust we consider, and in particular is inadequate for low densities in the outer layers, where one should account for degenerate electrons. In the current work we considered only a single polytropic EOS, adequate to describe degenerate neutrons, as this allowed us to obtain models for higher accreted masses. Polytropic models for degenerate electron EOSs (both relativistic and non-relativisitc) were also explored, but displayed a non-linear relation between the quadrupolar and spherical density perturbations already for values of the accreted mass much lower than those of interest, and did not allow to extrapolate the behaviour of the quadrupole to higher densities. Future work will aim to consistently model magnetic deformations and crustal reactions, with a realistic equation of state, in order to accurately estimate the induced quadrupole.

Finally we note that we have only considered deformations of the magnetic field and corresponding density profile. However the early evolution of the magnetic field due to Hall drift can lead to small scale structures that can persist for timescales greater than $10^{6}$ yrs, leading to sizable deformations Suvorov et al. (2016). If these structures can persist, even partially, in the crust of accreting neutron stars, they will lead to additional deformations at high density and deform deep reaction layers, possibly allowing for larger mountains and stronger gravitational wave emission.

Acknowledgements

We acknowledge support from the Polish National Science Centre grant SONATA BIS 2015/18/E/ST9/00577. Partial support comes from PHAROS, COST Action CA16214. TB and NS acknowledges support of the FNP grant TEAM/2016-3/19. BH wishes to thank J. L. Zdunik for interesting discussions and for sharing a preliminary version of the data in table (1). The authors are very grateful to D.I. Jones for useful comments and suggestions. NS wishes to thank N.  Andersson for insightful discussion and comments, and D. Bhattacharya for guidance and assistance in the initial development of the code. This document has been assigned LIGO document number P1900232