The physics beyond the standard model (BSM) has not appeared in extensive experimental and observational searches. Although the problems in the standard model strongly suggest the existence of a more fundamental theory, its detail is still unclear. In such a situation, it becomes more and more important to further extend the area of new physics search beyond the conventional ones. In this dissertation, we study the thermal evolution of neutron stars (NSs) as a probe of new physics.The NS cooling is well studied and the current standard theory is successful in explaining many of the observed surface temperatures of NSs. If a new particle is light enough to be created in NSs, and has small interactions with ordinary matter, its emission can be an extra cooling source. If this enhancement of cooling spoils the success of the standard cooling, the model should be excluded. An important candidate of such a new particle is axion, which dynamically solves the strong CP problem of the standard model. The axion-nucleon coupling induces the axion emission from the NS core, and it a ects the cooling of young NSs. The stringent constraint is obtained by comparing the predicted cooling curve to the observed surface temperature of the NS in the supernova remnant of Cassiopeia A (Cas A). We constrain the axion models by taking account of the temperature evolution in the whole life of the Cas A NS. The resultant limit on the axion decay constant is ≳ 10 8 GeV, which is comparable to the existing limit from SN1987A.The NS heating provides yet another example of probing new physics in NSs. In particular, it is known that the dark matter (DM) accretion leads to the heating of old NSs. It occurs through the DM accumulation in NSs by the scattering with nucleons, and their subsequent annihilation inside the star. The surface temperatures of old NSs ( ≳ 10 7 yr) are predicted to be = (2 − 3) × 10 3 K, which is stark contrast to the standard cooling theory. Thus the measurement of the surface temperatures of old NSs can provide a hint/constraint for DM models. This conclusion is, however, drawn with the assumption that there is no other heating source. In fact, it is known that the non-equilibrium beta process induces the late-time heating through the imbalance of chemical potentials among nucleons and leptons: this is called rotochemical heating and is inevitable for pulsars. The rotochemical heating typically predicts ∼ 10 5−6 K for old NSs, and hence can be much stronger than the DM heating. This possibility has been overlooked in the studies of DM heating. In this dissertation, we address the condition in which DM heating dominates the rotochemical heating. For that purpose, we rst perform detailed analysis of the rotochemical heating; extending the previous studies, we perform numerical calculation with both the proton and neutron pairing gaps, and compare the predictions to the observations. Then we include the DM heating, and investigate whether the DM heating is visible or not, varying the NS parameters. We nd that the DM heati...