An immersed boundary method for the thermo–fluid–structure interaction in rarefied gas flows is presented. In this method, the slip model is incorporated with the penalty feedback immersed boundary method to address the velocity and temperature jump conditions at the fluid–structure interface in rarefied gas flows within the slip-flow regime. In addition, the compressible flows governed by the Navier–Stokes equations are solved by using a high-order finite difference method; the elastic solid is solved by using the finite element method; the fluid and solid dynamics are solved independently, and the thermo–fluid–structure interaction is achieved by using a penalty feedback method in a partitioned way. To model the local rarefaction in the supersonic flow, an artificial viscosity is proposed by introducing the local Knudsen number to diffuse the sharp transition at the shock wave front. Several validations are conducted: the Poiseuille flow in a channel, the flow around a two-dimensional airfoil, a moving square cylinder in a channel, the flow around a sphere, and a moving sphere in quiescent flow. The numerical results from the present method show very good agreements with the previous published data obtained by other methods, confirming the good ability of the proposed method in handling the thermo–fluid–structure interaction in both weakly and highly compressible rarefied gas flows. Inspired by the micro/unmanned aerial vehicles in Martian exploration, the proposed method is applied to the aerodynamics of a flapping wing in rarefied gas flows in both two-dimensional and three-dimensional spaces to demonstrate the versatility of the proposed method for modeling flows involving large deformation and fluid–structure interaction.