Population synthesis studies constitute a powerful method to reconstruct the birth distribution of periods and magnetic fields of the pulsar population. When this method is applied to populations in different wavelengths, it can break the degeneracy
in the inferred properties of initial distributions that arises from single-band studies. In this context, we extend previous works to include $X$-ray thermal emitting pulsars within the same evolutionary model as radio-pulsars. We find that the cumulative distribution of the number of X-ray pulsars can be well reproduced by several models that, simultaneously, reproduce the characteristics of the radio-pulsar distribution. However, even considering the most favourable magneto-thermal evolution models with fast field decay, log-normal distributions of the initial magnetic field over-predict the number of visible sources with periods longer than 12 s. We then show that the problem can be solved with different distributions of magnetic field, such as a truncated log-normal distribution, or a binormal distribution with two distinct populations. We use the observational lack of isolated NSs with spin periods P>12 s to establish an upper limit to the fraction of magnetars born with B > 10^{15} G (less than 1%). As future detections keep increasing the magnetar and high-B pulsar statistics, our approach can be used to establish a severe constraint on the maximum magnetic field at birth of NSs.
