# Thermal conductivity and impurity scattering in the accreting neutron star crust

###### Abstract

We calculate the thermal conductivity of electrons for the strongly correlated multi–component ion plasma expected in the outer layers of neutron star’s crust employing a Path Integral Monte Carlo (PIMC) approach. This allows us to isolate the low energy response of the ions and use it to calculate the electron scattering rate and the electron thermal conductivity. We find that the scattering rate is enhanced by a factor 2–4 compared to earlier calculations based on the simpler electron–impurity scattering formalism. This finding impacts the interpretation of thermal relaxation observed in transiently accreting neutron stars, and has implications for the composition and nuclear reactions in the crust that occur during accretion.

^{†}

^{†}preprint: INT-PUB-16-003

## I Introduction

Observations of transient cooling in accreting neutrons stars and magnetars Eichler and Cheng (1989); Rutledge et al. (2002) after outbursts has motivated recent work to model the thermal evolution of the outermost regions of the neutron star called the crust Shternin et al. (2007); Brown and Cumming (2009); Page and Reddy (2012, 2013). In these models the temporal structure of the x-ray light curves is set by the thermal conductivity of the crust. Here, the relevant density and the expected temperature is in the range of K. For densities , electrons are relativistic and degenerate and they dominate transport phenomena. Under these conditions, nuclei are ionized and form a crystalline solid, and electron conduction is limited by electron-ion scattering. Since the typical electron wavelength is comparable to the distance between ions, interference between electron scattering off different ions is important. Accounting for such interference under arbitrary ambient conditions in a multi-component plasma (MCP) is a challenging many-body problem because ions have large atomic number (), and their dynamics is strongly correlated by Coulomb interactions at low temperature. In this article we present the first quantum calculation of electron-ion scattering in MCPs and electronic transport properties. The results we present apply directly to the outer crust, however the technique we propose here and some insights also apply to matter at higher density in the inner crust where some neutrons drip out off nuclei to occupy states in the continuum Baym et al. (1971).

In its simplest form, crustal matter is a one component plasma (OCP) of ions with charge in the range . The Coulomb interaction between ions is weakly screened by degenerate electron gas, and is given by where is the fine structure constant, is the electron screening length, and are the electron Fermi momentum and Fermi velocity, respectively. We note that here and throughout this article we use natural units and set . The ground state of a one component plasma (OCP) is a BCC solid at low temperature, and electron scattering in this phase, including the effects of the single and multi-phonon processes, is well studied in terrestrial metals Ziman (1960). Calculations of the electron thermal conductivity of an OCP over the full range of temperatures of interest to neutron star astrophysics are now available Flowers and Itoh (1976); Baiko et al. (1998); Potekhin et al. (1999); Abbar et al. (2015). However, when several species of ions are present these methods cannot be applied directly.

In accreting neutron stars a diverse mix of nuclei is produced through rapid proton capture (rp-process) reactions at the surface Schatz et al. (2001). We can expect several nuclear species to continue to coexist deeper in the crust because reaction pathways needed to process the rp-process ashes to the ground state nucleus are blocked by large Coulomb barriers. Their evolution through electron capture reactions at shallow depths and pycno-nuclear fusion reactions deeper in the crust produces a complex multi-component mixture of ions with a wide range of Gupta et al. (2007, 2008); Steiner (2012). The description of electron scattering in such a multi-component plasma (MCP) at low temperature is the major goal of this study. To date electron scattering in MCP has only been studies in the classical limit where the De Broglie wavelength of the ions is small compared to the average inter-ion distance. In these earlier studies molecular dynamics was used to provide a quantitive description of the thermal conductivity in the high temperature limit when Horowitz et al. (2009); Daligault and Gupta (2009) where

(1) |

is the average plasma frequency. Here, for each nuclear species labelled by the subscript , is the number density, is the mass, is the charge, and the abundance where is the mean ion density.

At lower temperatures of relevance to neutron stars, quantum effects cannot be neglected a priori. To include them we use the Path Integral Monte Carlo (PIMC) method (for a review see Ceperley (1995) ) to obtain the ion-ion correlation functions needed to calculate electron scattering rates and present first results of quantum calculations of the thermal conductivity. For typical MCPs encountered in neutron stars, we find significant reduction of the thermal conductivity at low temperature compared to those obtained in earlier work based on treating the MCP as an OCP plus uncorrelated impurities Flowers and Itoh (1976); Baiko et al. (1998); Potekhin et al. (1999).

The article is organized as follows. In section II we review well-known results for the electron thermal conductivity and its relation to the ion-ion correlation function. We discuss the quantum calculation of this correlation function in Euclidean time using PIMC in section III and present our results in section V. Finally we summarize and conclude in section VI.

## Ii Thermal conductivity

In the crust where electrons are degenerate and weakly coupled, their thermal conductivity

(2) |

where is the heat capacity of the relativistic electron gas, is the electron mean free path, and and are the electron Fermi velocity and Fermi energy, respectively Ziman (1960); Flowers and Itoh (1976). The second equality is obtained by noting that the specific heat of a degenerate electron gas is where is the electron density and electron collision rate . Under typical conditions electron-electron collisions are negligible and the total scattering rate and in the following we will only consider the electron-ion scattering process. This rate can be written as Nandkumar and Pethick (1984)

(3) |

where

(4) |

is the characteristic collision frequency and

(5) |

is called the Coulomb logarithm ^{1}^{1}1Although in the plasma physics literature the Coulomb logarithm is usually defined without in the integrand, in the context of dense astrophysical plasmas, the quantity defined by Eq. 5 is also called the Coulomb logarithm.. Here

(6) |

and , is the electron Fermi–momentum, is the Thomas–Fermi wave–vector and

(7) |

is the structure factor for the thermal conductivity Nandkumar and Pethick (1984) which contains all the information about ion–ion correlations. is the dynamic structure factor with contributions from elastic Bragg scattering removed because this contributes to the electron ground state wave-function and leads to their band structure but does not contribute to transport properties. denotes the average over the direction of unit vector and the function

(8) |

incorporates the final state blocking of electrons and detailed balance that ensures typical energy transfers is of the order of the temperature. Finally, the second term in parenthesis incorporates energy exchanging small angle scattering contributions to the thermal conductivity.

The charge–charge dynamic structure factor is the Fourier transform of the correlation function

(9) |

where the factor ensures the correct normalization and denote thermal averages at a temperature . The charge density operator is

(10) |

with and the charge and position of the i ion at time .

At high temperatures the bulk of the response is expected in the region where . Here and . However when this approximation fails and dynamical information is necessary to calculate . With decreasing temperature, electron scattering only probes the response at small of order the temperature and it is imperative to identify the strength of at to calculate the thermal conductivity.

In the astrophysical context the MCP is often approximated as a perfect crystal with ions of charge at the lattice sites plus a randomly distributed impurity charge Flowers and Itoh (1976); Itoh and Kohyama (1993). In this picture the total collision rate has two separate contributions . The first contribution is due to the absorption or emission of phonons by electrons. These processes are inelastic and consequently suppressed at low temperature. In contrast, impurity scattering is elastic and the temperature independent scattering rate is given by

(11) |

where

(12) |

is called impurity–parameter and the Coulomb logarithm for uncorrelated impurities is

(13) | |||||

(14) |

In the neutron star context one finds that even for modest values of the , for typical temperatures in the range K in the denser regions of the crust where Brown and Cumming (2009).

The above mentioned approach to describe electron scattering in the MCP is only approximate. It neglects correlations between minority and majority species and correlations between minority species can also be important as nuclei with small can cluster Horowitz et al. (2009). In the vicinity of an impurity we expect static distortions of the majority lattice. This would induce a non-periodic component to the Coulomb field which is bigger than the field associated with effective impurity charge . To account for these effects we require a method that can capture the dynamics of all species of ions on equal footing. In the classical limit, the dynamic structure function of the MCP has been calculated using Molecular Dynamics (MD), while the static structure function of MCP has been calculated also with Classical Monte Carlo (CMC) methods Horowitz et al. (2009); Daligault and Gupta (2009); Abbar et al. (2015). In what follows we describe the PIMC technique needed to perform the quantum calculation of the response of a MCP.

## Iii Path Integral Monte Carlo simulations

The Path Integral Monte Carlo method is an exact many–body technique to calculate equilibrium properties of strongly interacting quantum particles with either Boltzmann or Bose statistics Ceperley (1995). Because we expect ions to be spatially localized by strong Coulomb interactions even at low temperature when their thermal De Broglie wavelength are large, its a good approximation to treat them as Boltzmann particles. There are two major advantages to using PIMC instead of CMC. First, it naturally incorporates zero-point motion of the ions and second, it is possible to obtain some dynamical information about the system by computing the euclidean (imaginary–time) correlation function

(15) |

where is the imaginary–time interval. The latter will turn out to be particularly useful in describing at . As discussed earlier the elastic response with provides the temperature independent contribution to electron scattering and dominates at low temperature.

The traces in Eq. (15) are performed over the –dimensional configuration space of the ions and the full path is split into slices. This procedure allows reliable approximations for the density matrices at the higher (inverse) temperature . For the conditions explored in this study, we found the primitive approximation very accurate by comparing it to the exact two–particle density matrix obtained in the Feynman–Kac approach (see Ceperley (1995) for details). The final result of the PIMC calculation is at discrete values of imaginary–time and for a large number of momentum transfers compatible with the periodic boundary conditions of the simulation box, ie. with the length of the box and integers.

The only systematic bias present in the current computation is due to finite–size effects which in this case are mostly caused by the long–range nature of the interaction. However, due to screening a detailed resolution of long distance effects through Ewald summations Ewald (1917) is unnecessary. We find that summing over image charges in two nearest neighbor cells in all directions is adequate (this corresponds to including a total cells surrounding the simulation box). We have checked the converge of this procedure for both energies and structure factors using systems composed of different number of ions ranging from to and found convergence when in good agreement with earlier findings in Ref. Abbar et al. (2015). In what follows we present results obtained with systems with .

## Iv Multi-component Plasmas in Accreting Neutron Stars

In accreting neutron stars rapid proton capture reactions can produce a diverse mix of nuclei at the surface with Schatz et al. (1999, 2001). As this mixture is incorporating into the outer crust phase separation and electron capture reactions purify the mix somewhat and the impurity parameter in the outer crust is expected to be in the range (cf. Gupta et al. (2007); Horowitz et al. (2007)). In this study we approximate the complicated composition found in Gupta et al. (2007) with a one component plasma with a charge and two- and three–component system labelled MCP1, MCP2 in Tab. 1 to mimic and . In addition we considered a two–component system, labelled MCP3, with motivated by findings in Ref. Brown and Cumming (2009); Turlione et al. (2015). The ion charges, their abundances denoted by and the impurity parameters are shown in Tab. 1. In all cases the mass density has been fixed to .

Z | A | x | ||

OCP | ||||

MCP1: | ||||

MCP2: | ||||

MCP3: | ||||

As can be expected at low temperature we find a large number of metastable states with large barriers and setting up the initial conditions is quite challenging. A careful analysis of the equilibrium configuration in a low temperature MCP is beyond the scope of this work, and would require an annealing algorithm to evolve to the low temperatures structure starting from a high temperature liquid configuration. Such calculations have been performed in the past using molecular dynamics simulations Horowitz et al. (2007, 2009); Horowitz and Berry (2009) where it was found that the ground state was a BCC crystal composed by ions with large charge at regular lattice sites and low–Z ions where found to occupy interstitial regions. This is not surprising since the Coulomb interaction is quite soft and allows for the diffusion in the solid phase that ensures amorphous structures to relax to crystalline state Hughto et al. (2011).

Motivated by these findings we initialize the ions on lattice–sites of a perfect BCC crystal which is then distorted by applying random displacement to all the particles. The system is allowed to relax to its ground–state after which statistics for observable start to be taken. In a multi–component system we also choose the type of ion on a given lattice at random, and during equilibration ions are allowed to interchange location with others of different species. As can be seen from the presence of BCC Bragg peaks in the structure factors in Fig. 1, the equilibrium configuration attained with this procedure is always an ordered BCC crystal phase. When ions with very small Z are present as in the MCP2, we found that in our simulations they remained as substitutional impurities on lattice sites. This may be an artifact of our simple initialization procedure and it would be interesting to properly explore annealing in future using replica–exchange methods Marinari and Parisi (1992); Sugita and Okamoto (1999). However, since this only affects the distribution of ions with very small species we do not expect their dynamics to dominate the charge-charge correlation function.

## V Results

We first present results for the static structure factor of the OCP and the MCPs for the compositions shown in Tab. 1 at three different temperatures MeV corresponding to . Calculations employing Classical Monte Carlo (where the path–integral is restricted to a single time–slice) were also performed to asses the importance of zero–point motion. For the OCP we confirm earlier results in Abbar et al. (2015) where it was found that the static structure function is underestimated in the classical case when .

In Fig. 1 we show for the OCP and MCP2 systems at the lowest temperature . It is apparent from the appearance of Bragg peaks that the underlying BCC structure persists in this MCP in agreement with earlier results obtained using MDHorowitz and Berry (2009).

For the OCP it is well known that the one phonon approximation provides a good description of S(q) at low temperature and that multi-phonon effects can be accounted for within the Harmonic Approximation. In Baiko et al. (1998), a simple analytic formula was developed based on the Harmonic Approximation (HA) to describe the static structure factor of an OCP and is shown as the solid black curve labelled OCP-FIT in the figure. Explicitly this is given by

(16) |

where the Debye–Weller factor was found to be of the following simple analytic form

(17) |

to describe systems over a wide range of temperature Baiko and Yakovlev (1995). Here and are the moments of the phonon spectrum. For a BCC lattice and Pollock and Hansen (1973). For the MCP2 the simple impurity model predicts

(18) |

since the contribution due to uncorrelated random impurities is additive. This is shown as the black-dashed curve in the figure and is labelled (OCP+imp-FIT).

Apart from sharp features coming from the Bragg peaks, the analytical fit is on average in good agreement with our numerical result for . However, the comparison between our results for the MCP2 and those obtained using the OCP + impurity fit show important differences. Impurity effects included through Eq. 11 significantly underestimates the contribution to and this is in qualitative agreement with results found using MD simulations at higher temperatures Horowitz et al. (2009); Horowitz and Berry (2009).

We now turn to discuss how we obtain informations about the dynamical structure factor from the Euclidean response function . From Eq. (15) one can in principle try to perform a numerical Laplace transform inversion starting with the PIMC data to obtain the frequency dependent response. However, as is well known McWhirter and Pike (1978); Talbot (1979), this numerical inversion is an ill-posed problem in the sense that arbitrarily small perturbations in the initial data can give rise to arbitrarily large deviations in the final answer. We attempted this inversion with different techniques and regularization schemes (Magierski and Wlazlowski (2012); Roggero et al. (2013) and references therein) with little success. This is because is very sensitive to response in the low energy at low temperatures, whereas the response obtained from numerical inversion is expected to be accurate in the region where the spectrum has maximum strength, and these two regions do not overlap significantly for most momentum transfers.

Due to the above mentioned problems, we choose instead to constrain models for the dynamical structure factor with PIMC data for . First, for the OCP we calculate the Euclidean response function using obtained in the HA and taking the Laplace transform as defined in Eq. (15). In Fig. 2 we compare the (normalized) euclidean response for charge–charge fluctuations obtained in this way with the predictions of PIMC for a representative fixed momentum transfer at three different temperatures for the OCP system. The excellent agreement between PIMC results and HA model is expected and confirms that the harmonic approximation works quantitatively. The striking feature is the rapid decrease in for large at low temperature. This simply reflects the fact that in the OCP the low energy response due to phonon excitations is highly suppressed at low temperature and demonstrates that PIMC data for at provides powerful means to extract the response at very low energy when .

The Euclidean correlation function in the MCP is shown in Fig. 3 for a momentum transfer at . A comparison between the results for the OCP (shown as black squares) and the MCP at shows as expected that imperfections in the lattice produce a significant response at low energy. The presence of excitations at can be directly inferred by the magnitude of the correlator at the largest imaginary–time separation : in the OCP, goes to zero at large imaginary–time indicating the absence of low energy excitations for , whereas in the MCP systems the persistence of correlations indicate cleanly the presence of significant strength at these low energies. This low energy response can be thought of as compositional modes with long relaxation times, corresponding to nearly static distortions of the periodic potential expected in the OCP. Molecular dynamics studies by Caballero et al. Caballero et al. (2006) also found additional strength at very low energy in MCP.

The absence of structure in for seen in Fig. 3 indicates that there is a clear separation between the low and high energy responses. This motivates us to write the dynamical structure factor at low temperature as

(19) |

where we separate the contributions at finite frequency due to excitation of phonons and an elastic contribution due to nearly static lattice imperfections generated by the presence of impurities. Taking its Laplace transform we obtain

(20) |

In general the extraction of will rely on a model for that can capture the dependence of the full Euclidean correlator, however since in the low temperature limit

(21) |

This identification is one important result of our study and provides a simple and robust strategy to calculate the thermal conductivity of complex mixtures encountered in accreting neutron stars for . In practice we find that for the contribution of the phonons is small and there is no difference between using Eq. (21) and Eq. (20) to extract independently of the model for . We tried calculating using for the average OCP and for the linear mixing model Daligault and Gupta (2009); Potekhin et al. (1999) for the MCP in the harmonic approximation and found negligible differences. In the following we will denote results obtained using Eq. (20) or Eq. (21) as "EUC".

For high temperature when we have already noted that and dynamical information contained in is not necessary to calculate the thermal conductivity. At moderate temperature the situation is complicated because the elastic contribution from phonons can be significant for large and the extraction of the contribution due to impurities will have some model dependence. Calculations of in MCP using molecular dynamics in Ref. Caballero et al. (2006) showed that the finite frequency contribution to the dynamical structure factor is very similar to the expected in the OCP with . This allows us to approximate the impurity contribution at moderate temperatures as

(22) |

where is calculated using PIMC. We note that Eq. (22) is remarkably accurate as its predictions at low temperature are consistent with the model–independent extraction obtained using Eq. (21).

In Fig. 4 we plot the Coulomb logarithm associated with the impurity contribution and defined by

(23) |

for various temperatures and for three multi-component mixtures presented in Tab. 1.

It is interesting to note that is approximately linear in and this allows us to define an effective impurity parameter

(24) |

The value extracted for is found to be temperature dependent: we obtain for and asymptotes to in the low temperature limit. To explore the temperature dependence we plot as a function of the Coulomb coupling parameter in Fig. 5. The figure also shows results obtained from molecular dynamics simulations of mixtures with Horowitz et al. (2009) (black diamonds in the figure) and Horowitz and Berry (2009) (orange diamonds) performed at a much higher density g/cm. It is important to stress that, despite the higher density, these calculations where targeted to describe the outer crust, and therefore neglect effects of finite size of the ions as well as the presence of interstitial neutrons. There is good agreement between our results and those obtained with molecular dynamics at smaller value of but the trend with increasing is different. The origin of this difference is unclear and warrants further work.

## Vi Summary and Conclusions

We have calculated the Euclidean charge-charge correlation of MCP using PIMC simulations and used it to extract the low energy response needed to calculate the electron thermal conductivity of high density matter at low temperatures encountered in neutron stars. The behavior of the Euclidean correlation function at large imaginary time was found to be dominated by slowly varying compositional modes associated with impurity induced lattice distortions. These contribute to the low energy strength in the ion dynamic structure and dominates electron scattering when . The role of impurity scattering in MCP has been studied in previous work in the high temperature classical regime using MD simulations Caballero et al. (2006). Our study is the first to extend the results to lower temperatures of relevance to neutron stars and includes quantum effects. The results we obtain are in agreement with the MD results at higher temperature which found that the simple impurity scattering formula underestimates the low energy response Horowitz et al. (2009); Horowitz and Berry (2009). However, we find that at low temperature this enhancement is substantially larger and implies that when the electron collision rate in typical MCP encountered in accreting neutron stars is about factor of 4 larger than earlier estimates based on the simple impurity scattering formula.

We have performed calculations for three MCP characterized by the impurity parameter and at various temperatures and find that the electron scattering rate increases linearly with . This allowed us to define a temperature dependent effective impurity parameter which can be used to calculate the electron thermal conductivity using the simple impurity formula in Eq. (11). In our simulations we found that the ions are highly localized and the Bragg peaks associated with a BCC lattice persists. This suggests that electron scattering is predominantly due to distortions of the lattice rather than due to scattering of randomly distributed impurity charge and is a likely explanation for the enhanced response at low energies that we observe. When some fraction of the impurities occupy interstitial spaces through thermal motion these distortions can be screened and offer an explanation for the reduction in the at higher temperature.

Finally, we comment on the implications of our study for the interpretation of observed thermal relaxation in accreting neutron stars. Models that best fit observations require a relatively large thermal conductivity in the inner crust with Shternin et al. (2007); Brown and Cumming (2009); Page and Reddy (2013); Turlione et al. (2015). Although our results do not directly apply to the inner crust because of the presence of dripped neutrons, it is reasonable to expect that qualitative trends will be similar since the Coulomb interaction will continue to be the dominant forces between ions. The enhanced scattering rates implied by indicates that the impurity parameter in the inner crust can be a factor of a few smaller than . With future observations and improved modeling it may be possible to obtain useful constraints on the low temperature conductivity of the outer crust. Our calculations combined with such input will provide useful constraints on the composition of the outer crust of an accreting neutron star.

###### Acknowledgements.

We would like to thank S. Abbar, J. Carlson, H. Duan and F. Pederiva for useful discussions. The work of S. R. was supported by the DOE Grant No. DE-FG02-00ER41132 and by the Joint Institute for Nuclear Astrophysics (JINA-CEE). The work of A. R. was supported by NSF Grant No. AST-1333607. Most of the intensive computations have been performed at NERSC thank to a Startup allocation.## References

- Eichler and Cheng (1989) D. Eichler and A. F. Cheng, Astrophys. J. 336, 360 (1989).
- Rutledge et al. (2002) R. E. Rutledge, L. Bildsten, E. F. Brown, G. G. Pavlov, V. E. Zavlin, and G. Ushomirsky, Astrophys.J. 580, 413 (2002), astro-ph/0108125 .
- Shternin et al. (2007) P. S. Shternin, D. G. Yakovlev, P. Haensel, and A. Y. Potekhin, Monthly Notices of the RAS 382, L43 (2007), arXiv:0708.0086 .
- Brown and Cumming (2009) E. F. Brown and A. Cumming, Astrophys. J. 698, 1020 (2009), arXiv:0901.3115 [astro-ph.SR] .
- Page and Reddy (2012) D. Page and S. Reddy, Thermal and transport properties of the neutron star inner crust, edited by C. Bertulani and J. Piekarewicz, Space Science, Exploration and Policies (NOVA, 2012) Chap. 14, arXiv:1201.5602 [nucl-th] .
- Page and Reddy (2013) D. Page and S. Reddy, Phys.Rev.Lett. 111, 241102 (2013).
- Baym et al. (1971) G. Baym, H. A. Bethe, and C. Pethick, Nucl. Phys. A175, 225 (1971).
- Ziman (1960) J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids. (Oxford University Press, 1960).
- Flowers and Itoh (1976) E. Flowers and N. Itoh, ApJ 206, 218 (1976).
- Baiko et al. (1998) D. A. Baiko, A. D. Kaminker, A. Y. Potekhin, and D. G. Yakovlev, Phys. Rev. Lett. 81, 5556 (1998).
- Potekhin et al. (1999) A. Y. Potekhin, D. A. Baiko, P. Haensel, and D. G. Yakovlev, Astron. Astrophys. 346, 345 (1999).
- Abbar et al. (2015) S. Abbar, J. Carlson, H. Duan, and S. Reddy, ArXiv e-prints (2015), arXiv:1503.01696 [astro-ph.HE] .
- Schatz et al. (2001) H. Schatz, A. Aprahamian, V. Barnard, L. Bildsten, A. Cumming, M. Ouellette, T. Rauscher, F.-K. Thielemann, and M. Wiescher, Phys. Rev. Lett. 86, 3471 (2001).
- Gupta et al. (2007) S. Gupta, E. F. Brown, H. Schatz, P. Moeller, and K.-L. Kratz, The Astrophysical Journal 662, 1188 (2007).
- Gupta et al. (2008) S. S. Gupta, T. Kawano, and P. Moller, Phys. Rev. Lett. 101, 231101 (2008), arXiv:0811.1791 [astro-ph] .
- Steiner (2012) A. W. Steiner, Phys. Rev. C85, 055804 (2012), arXiv:1202.3378 [nucl-th] .
- Horowitz et al. (2009) C. J. Horowitz, O. L. Caballero, and D. K. Berry, Phys. Rev. E 79, 026103 (2009).
- Daligault and Gupta (2009) J. Daligault and S. Gupta, The Astrophysical Journal 703, 994 (2009).
- Ceperley (1995) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- Nandkumar and Pethick (1984) R. Nandkumar and C. J. Pethick, Mon. Not. Roy. Astro. Soc. 209, 511 (1984).
- Itoh and Kohyama (1993) N. Itoh and Y. Kohyama, ApJ 404, 268 (1993).
- Ewald (1917) P. Ewald, Annal. Phys. 54, 557 (1917).
- Schatz et al. (1999) H. Schatz, L. Bildsten, A. Cumming, and M. Wiescher, Astrophys. J. 524, 1014 (1999), arXiv:astro-ph/9905274 [astro-ph] .
- Horowitz et al. (2007) C. J. Horowitz, D. K. Berry, and E. F. Brown, Phys. Rev. E 75, 066101 (2007).
- Turlione et al. (2015) A. Turlione, D. N. Aguilera, and J. A. Pons, A&A 577, A5 (2015).
- Horowitz and Berry (2009) C. J. Horowitz and D. K. Berry, Phys. Rev. C 79, 065803 (2009).
- Hughto et al. (2011) J. Hughto, A. S. Schneider, C. J. Horowitz, and D. K. Berry, Phys. Rev. E 84, 016401 (2011).
- Marinari and Parisi (1992) E. Marinari and G. Parisi, EPL (Europhysics Letters) 19, 451 (1992).
- Sugita and Okamoto (1999) Y. Sugita and Y. Okamoto, Chemical Physics Letters 314, 141 (1999).
- Baiko and Yakovlev (1995) D. A. Baiko and D. G. Yakovlev, Astron. Lett. 21, 702 (1995).
- Pollock and Hansen (1973) E. L. Pollock and J. P. Hansen, Phys. Rev. A 8, 3110 (1973).
- McWhirter and Pike (1978) J. McWhirter and E. Pike, J. Phys. A 11(9), 1729 (1978).
- Talbot (1979) A. Talbot, J. Inst. Maths. Applics. 23(1), 97 (1979).
- Magierski and Wlazlowski (2012) P. Magierski and G. Wlazlowski, Comput. Phys. Commun. 183, 2264 (2012).
- Roggero et al. (2013) A. Roggero, F. Pederiva, and G. Orlandini, Phys. Rev. B 88, 094302 (2013).
- Caballero et al. (2006) O. L. Caballero, C. J. Horowitz, and D. K. Berry, Phys. Rev. C 74, 065801 (2006).