We study tunneling spectroscopy of subgap Andreev states in a superconducting system and discuss the general situation when the discrete nature of these levels is relevant and thus the standard semiclassical result for tunneling conductance being proportional to the density of states is not applicable. If the tunneling coupling is weak, individual levels are resolved and the conductance G(V ) at low temperatures is composed of a set of resonant Lorentz peaks which cannot be described within perturbation theory over tunneling strength. We establish a general formula for the peak widths and heights and show that the width of any peak scales as normal-state tunnel conductance, while its height is 2e2?h?1 and depends only on contact geometry and the spatial profile of the resonant Andreev level.

We also establish an exact formula for the single-channel conductance that takes the whole Andreev spectrum into account, and use it to study the interference of Andreev reflection processes through different levels. We study tunneling conductance at finite bias G(eV >0) for a system with a pair of almost decoupled Majorana fermions and derive the conditions for the 'universal' zero-bias peak with the height 2e2/h to be observed in a realistic system which always hosts an even number of Majorana fermions.