A novel methodology, referred to as Auxiliary Domain Method (ADM), allowing for a very efficient solution of nonlinear reliability problems is presented. The target nonlinear failure domain is first populated by samples generated with the help of a Markov Chain. Based on these samples an auxiliary failure domain (AFD), corresponding to an auxiliary reliability problem, is introduced. The criteria for selecting the AFD are discussed. The emphasis in this paper is on the selection of the auxiliary linear failure domain in the case where the original nonlinear reliability problem involves multiple objectives rather than a single objective. Each reliability objective is assumed to correspond to a particular response quantity not exceeding a corresponding threshold. Once the AFD has been specified the method proceeds with a modified subset simulation procedure where the first step involves the direct simulation of samples in the AFD, rather than standard Monte Carlo simulation as required in standard subset simulation. While the method is applicable to general nonlinear reliability problems herein the focus is on the calculation of the probability of failure of nonlinear dynamical systems subjected to Gaussian random excitations. The method is demonstrated through such a numerical example involving two reliability objectives and a very large number of random variables. It is found that ADM is very efficient and offers drastic improvements over standard subset simulation, especially when one deals with low probability failure events.