  Obliquity variability could play an important role in the climate and
habitability of a planet. Orbital modulations caused by planetary companions
and the planet's spin axis precession due to the torque from the host star may
lead to resonant interactions and cause large-amplitude obliquity variability.
Here we consider the spin axis dynamics of Kepler-62f and Kepler-186f, both of
which reside in the habitable zone around their host stars. Using
{\emph{N}}-body simulations and secular numerical integrations, we describe
their obliquity evolution for particular realizations of the planetary systems.
We then use a generalized analytic framework to characterize regions in
parameter space where the obliquity is variable with large amplitude. We find
that the locations of variability are fine-tuned over the planetary properties
and system architecture in the lower-obliquity regimes ($\lesssim 40^\circ$).
As an example, assuming a rotation period of 24 hr, the obliquities of both
Kepler-62f and Kepler-186f are stable below $\sim 40^\circ$, whereas the
high-obliquity regions ($60^\circ - 90^\circ$) allow moderate variabilities.
However, for some other rotation periods of Kepler-62f or Kepler-186f, the
lower-obliquity regions could become more variable owing to resonant
interactions. Even small deviations from coplanarity (e.g. mutual inclinations
$\sim 3^\circ$) could stir peak-to-peak obliquity variations up to $\sim
20^\circ$. Undetected planetary companions and/or the existence of a satellite
could also destabilize the low-obliquity regions. In all cases, the
high-obliquity region allows for moderate variations, and all obliquities
corresponding to retrograde motion (i.e. $> 90^\circ$) are stable.
