The local reconstruction of conservative fluxes from finite element solutions plays a crucial role in numerical analysis, particularly in applications such as a posteriori error estimation and the enforcement of local conservation properties in continuum mechanics. These topics have been extensively investigated in the literature for problems discretized on fitted meshes. In this work, we focus on extending flux reconstruction techniques to the unfitted mesh setting. More precisely, we consider the CutFEM method, which allows the numerical treatment of interfaces that are not aligned with the background mesh. Within this framework, we propose a local reconstruction of the flux in an immersed Raviart–Thomas space. and perform a flux-based a posteriori error analysis together with the corresponding adaptive mesh refinement strategy. We study a second-order elliptic interface problem in which the interface is not aligned with the mesh, with standard transmission conditions. The diffusion coefficients are assumed to be discontinuous across the interface. The numerical approximation is based on continuous, piecewise linear finite elements, and the interface is handled using the CutFEM method [2]. The flux reconstruction follows the approach developed in [1] for the Poisson problem on fitted meshes. The idea is to introduce a hybrid mixed formulation, equivalent to the discrete primal formulation. This formulation involves a Lagrange multiplier defined on the mesh edges, which can be computed locally and subsequently used to correct the numerical flux (see [3]). We next define a conservative flux in the lowest-order immersed Raviart–Thomas space of [5], whose elements satisfy the transmission condition across the interface in a strong sense, while also weakly enforcing continuity of the solution. Finally, we introduce a new flux-based a posteriori error estimator that contains both volume and interface contributions. We establish sharp reliability and local efficiency results [4], with particular attention paid to the dependence of the efficiency constants on the diffusion coefficients and on the mesh / interface configuration. Numerical experiments are presented, illustrating the robustness of the approach.