What is the core temperature of neutron star as a function of time? I wish to know about the temperature of the bulk (not just the surface) of a neutron star. I am especially interested in the highest temperature part which I assume is the core, but would want to know how much of the star that is. I mainly want to know about the early part (thousands of years) after formation, but any information would be helpful.
I have read statements about billions of kelvin and then millions of kelvin, about MeV and then keV, but I have not managed to get clarity about timescales. A paper from 1972 said the temperature might remain $\sim 10^9$ K for 'a thousand' years. Some recent websites say it cools to $10^6$ in 'a few years'. Which is right? If the earlier statement is wrong is it because neutrino processes were not understood correctly back in the 1970s?
Also you often see statements that neutron stars have this temperature or that temperature, but surely if the star is cooling then what temperature is being referred to? When in the cooling are we and how quickly is it falling? There is a graph of surface temperature vs time in Fig 1 of
https://arxiv.org/abs/astro-ph/0402143
for stars of not too high mass which cool more slowly. The plotted line starts at $10^{6.4}$ K. In the same paper after Fig 3 it says "The internal temperature drops to $T_i \sim 1.5 \times 10^8$ K in $t = 10^5$ yrs". So this suggests the ratio of internal to surface temperature can be about 200. Can that value be assumed at earlier times when the star is hotter? (Statements about surface temperature are not much use if you want to know what is going on inside.)
What is a good ballpark figure for the ratio between core and surface temperature?
 A: It's a very broad and interesting question. I am only going to provide a summary here and refer to some literature (since I know you are an "expert user") and since there aren't any really simple or quick answers to your questions.
A number of important points to consider:
Neutron stars are born very hot - probably of order $10^{11}$ K and there is some evidence for that from the observation of neutrinos from SN 1987A.
At these temperatures, the neutrons, and protons in the core are partially degenerate and the core varies from being opaque to neutrinos for a few seconds, to partially opaque and then effectively transparent as the neutron star cools. The reason for this is that once the neutron star cools and becomes fully degenerate at $T < 10^{10}$ K, then neutrino capture or scattering is highly suppressed by the lack of available momentum states for the neutrons or protons. The opacity or otherwise of the core is crucial to the thermal history because it is only escaping neutrinos that allow the neutron star to cool - you can forget about photon cooling from the surface during the first 1000 years, it is totally ineffective by comparison because of (i) the small surface area and (ii) there is indeed a substantial temperature gradient in the surface regions which means that the "photospheric" temperature is about a factor of 100 lower than the interior temperature.
Heat transport within a neutron star can be very efficient. High levels of degeneracy mean the thermal conductivity is extremely high. Neutron star interiors should be more-or-less isothermal, other than in the very early stages of cooling when the various regions can cool at different (very high) rates. This is because of different levels of degeneracy and the density gradient in the star and cooling can be fast enough in the interior that thermal conductivity can't eliminate the temperature gradients fast enough to achieve isothermality.
The main cooling processes are the URCA and modified URCA processes. These are cycle of beta- and inverse-beta decay that result in the emission of (anti) neutrinos. The former and more basic process is vastly more efficient and could cool the neutron star interior in seconds if it could operate unimpeded. However, once degeneracy sets in, the URCA process is blocked by a lack of available momentum states and bystander particles are needed to conserve energy and momentum simultaneously - this is the modified URCA process.
Because this requires all participating particles to be within $\sim k_B T$ of their respective Fermi surfaces, it is relatively inefficient. Both processes have a very high temperature dependence ($T^6$ and $T^8$ respectively). There is also neutron-neutron bremsstrahlung, releasing neutrinos, which is less efficient again, but which may become important if even modified URCA is suppressed by superfluidity at lower temperatures (see below). However, the direct URCA process can continue to operate even in a highly degenerate gas as long as the n/p ratio is low enough. Since this ratio decreases with density, it means there is a threshold density inside the star, interior to which the direct URCA process and fast cooling can continue, but slower cooling occurs in regions outside this. Since interior densities are likely mass dependent, this introduces a mass-dependence to the cooling.
The heat capacity of a neutron star is very low once degeneracy is achieved within the first few seconds. Basically, only neutrons within $k_B T$ of their Fermi surface can be cooled. This, along with the ease with which neutrinos can escape, is the reason neutron stars cool down extremely rapidly. There is a twist though to do with superfluidity - this is a possibility once the interior cools below a critical temperature (probably around $k_B T_c \sim 1$ MeV) when neutrons can pair and this creates a sudden drop in the heat capacity at those temperatures. On the other hand it also suppresses the cooling neutrino cooling processes when $T\ll T_C$  because the pairs have to be broken before they can participate in the modified URCA process but enhances the neutrino cooling at $T \sim T_c$ because the neutron pairs can be formed and broken producing additional thermal neutrinos. This is a complication which isn't fully resolved in a theoretical sense.
Another complication is what is going on right at the core of the neutron star. If densities are high enough then it is possible that additional hadronic phases could emerge (things like heavy $\Lambda$ particles). These could decrease the n/p ratio to an extent that allows the direct URCA process to operate for longer than it otherwise would. There is also the possibility of quark matter if "asymptotic freedom" is achieved. Quark matter is also able to cool by the direct URCA process and so neutron stars featuring quark matter should cool (much) more rapidly. Since the central densities of neutron stars depend on their mass (in a way that depends on an uncertain equation of state) then it could be that neutron stars of different masses cool very differently.
Finally, you ask about the difference between the interior temperature and "surface temperature". The latter is what one attempts to observe with, for example, X-ray astronomy satellites, since the interior temperature is not directly observable at present. The ratio between the interior and surface temperature is something like a factor of 100-1000 and it is higher when the interior is hotter. It also depends to some extent on things like magnetic fields that can change the mechanical equilibrium in the outer layers and mess with heat transport mechanisms. A rough formula that can be used is
$$
\frac{T_{\rm eff}}{T_{\rm int}} = 9.3\times10^{-3} \left(\frac{T_{\rm int}}{10^{8} {\rm K}}\right)^{-0.45}\, , $$
where $T_{\rm eff}$ is the effective surface temperature and $T_{\rm int}$ is the interior temperature (from Shapiro & Teukolsky 1983).
So the above is a basic summary of the physics involved - as a (bit more advanced than) primer to the subject I highly recommend the authoritative review by Yakovlev & Pethick (2004) from which the following plots and timeline are taken.
Neutron stars are born with cores at $\sim 10^{11}$ K and spend their first minute or so as proto-neutron stars, opaque to neutrinos. Then for 10-100 years rapid neutrino cooling via the direct URCA process can actually make the core cooler than the crust - so while the core cools rapidly, the crust and surface temperatures cool more slowly. From 100-$10^5$ years, the interior becomes isothermal and neutrino cooling dominates over photon cooling. The details of how the interior (and hence surface) temperature behaves then depends greatly on the composition of the core, hence the mass of the neutron star, and the extent to which superfluidity operates in the neutron fluid. The plots below are typical, showing the evolution of surface temperature with time (right panel), while the left-hand panel is showing the efficiency normalisation for the cooling rate as a function of density in the neutron star interior.
As an orientation you can see from the equation above that when the surface temperature is $10^{5.5}$ K, that $T_{\rm int} \sim 10^7$ K, if $T_{\rm eff} \sim 10^6$ K then $T_{\rm int} \sim 10^8$ K and $T_{\rm eff} = 10^{6.5}$ K would correspond to $T_{\rm int} \sim 10^{9}$ K.
The upper curve is appropriate for low mass ($\sim 1.35M_{\odot}$) neutron stars made of entirely nucleonic matter and with no superfluidity and suppressed direct URCA cooling. The lower curves show what happens for higher mass, higher density neutron stars that are able to cool more directly. You can see that a wide range of interior temperatures at 1000 years is possible - anything from $\sim 2 \times 10^8$ K down to maybe $10^7$ K

