AbstractParticle swarm optimization, one of the modern global optimization methods, is attracting widespread interest because it overcomes the difficulties of conventional inversion techniques, such as trapping at a local minimum and/or initial model dependence. The main characteristic of particle swarm optimization is the large search space of parameters, which in a sense allows the exploration of the entire objective function space if the input parameters are properly chosen. However, in the case of a high‐dimensional model space, the numerical instability of the solution may increase and lead to unrealistic models and misinterpretations due to the sampling problem of particle swarm optimization. Therefore, smoothness‐constrained regularization techniques used for the objective function or model reduction techniques are required to stabilize the solution. However, weighting and combining objective function terms is partly a subjective process, as the regularization parameter is generally chosen based on some kind of criteria of how the smoothing constraints affect the data misfits. This means that it cannot be completely predefined but needs to be adjusted during the inversion process, which begins with the response of an initial model. In this paper, a new modelling approach is proposed to obtain a smoothness‐constrained model from magnetotelluric data utilizing multi‐objective particle swarm optimization based on the Pareto optimality approach without using a regularization parameter and combining several objective function terms. The presented approach was verified on synthetic models and an application with field data set from the Çanakkale–Tuzla geothermal field in Turkey. Findings from these analyses confirm the usefulness of the method as a new approach for all constrained inversions of geophysical data without the need to combine the objective function terms weighted by a regularization parameter.