Global stability analysis of the WASP-148 planetary system. We fixed all orbital parameters of the solution shown in Table 6 and varied the only two unconstrained parameters, both from the outer planet WASP148c: the longitude of its ascending node, Ω_c, and its inclination, i_c. The upper panel shows the whole (i_c, Ω_c) domain, whereas the lower panel zooms ino the most stable regions. The step size was 0.2° in the node and 0.5° in the inclination. For each initial condition, the system was integrated over 50 kyr and a stability criterion was derived with the frequency analysis of the mean longitude. Because of the dynamical invariant, the figure is symmetric with respect to the (i_c = 90°, Ω_c = 0) center. White dashed curves give the isolines of constant mutual inclination I = 10°, 20°, 30°, and 40°. The color scale corresponds to values between − 9 (black) and − 3 (red) for the decimal logarithm of the stability index D used in Correia et al. (2010). The red zones correspond to highly unstable orbits, while the dark blue region can be assumed to be stable on a billion-year timescale. The reduced-χ² level curves of the Newtonian fit are also plotted in white.