Let $f:\mathbb{R}^2\to\mathbb{R}$. Can we do Taylor approximation for only one variable $$f(x,y) \approx f(x_0,y) + \frac{\partial }{\partial x}f(x_0,y)(x-x_0) + \frac{1}{2}\frac{\partial^2}{\partial x^2}f(x_0,y)(x-x_0)^2$$
I thought that we can always do this, by treating the other variable $y$ as fixed, but I have some hesitations.
You can certainly form this expansion---after all, it's simply producing the Taylor polynomial in $x$ for each single-variable function $f(x, y_0)$---and then one can apply the usual results of single-variable Taylor series (error estimates, etc.).
Whether you want to do this depends on context, though; in particular, NB that the approximation you've given is not polynomial in $x,y$, and in some sense the purpose of the Taylor approximations is that they are polynomial, and hence often easier to work with than the original function.