  We test how metallicity variation (a background gradient and fluctuations)
affects the physics of local thermal instability using analytical calculations
and idealized, high-resolution 1D hydrodynamic simulations. Although the
cooling function ($\Lambda[T,Z]$) and the cooling time ($t_{\rm cool}$) depend
on gas temperature and metallicity, we find that the growth rate of thermal
instability is explicitly dependent only on the derivative of the cooling
function relative to temperature ($\partial \ln \Lambda/\partial \ln T$) and
not on the metallicity derivative ($\partial \ln \Lambda/ \partial \ln Z$). For
most of $10^4~{\rm K} \lesssim T \lesssim 10^7~{\rm K}$, both the isobaric and
isochoric modes (occurring at scales smaller and larger than the sonic length
covered in a cooling time [$c_s t_{\rm cool}$], respectively) grow linearly,
and at higher temperatures ($\gtrsim 10^7~{\rm K}$) the isochoric modes are
stable. We show that even the nonlinear evolution depends on whether the
isochoric modes are linearly stable or unstable. For the stable isochoric
modes, we observe the growth of small-scale isobaric modes but this is distinct
from the nonlinear fragmentation of a dense cooling region. For unstable
isochoric perturbations we do not observe large density perturbations at small
scales. While very small clouds ($\sim {\rm min}[c_st_{\rm cool}]$) form in the
transient state of nonlinear evolution of the stable isochoric thermal
instability, most of them merge eventually.
