Variation of the Lagrangian and the Noether current In Schwartz’s book, QFT and Standard Model, section 8.3.1, he writes

if we then let $\alpha$ be a function of $x$, the transformed $\mathcal L_0$ can only depend on $\partial_\mu \alpha$. Thus, for infinitesimal $\alpha(x)$,
$$\delta \mathcal{L}_{0}=\left(\partial_{\mu} \alpha\right) J_{\mu}+\mathcal{O}\left(\alpha^{2}\right)\tag{8.58}$$

for the Noether current defined as
$$J_{\mu}=\sum_{n} \frac{\partial \mathcal{L}}{\partial\left(\partial_{\mu} \phi_{n}\right)} \frac{\delta \phi_{n}}{\delta \alpha},\tag{8.57}$$
and that

in scalar QED, with $A_{\mu}=0, \mathcal{L}_{0}=\left(\partial_{\mu} \phi\right)^{\star}\left(\partial_{\mu} \phi\right)-m^{2} \phi^{\star} \phi$ and
$$\delta \mathcal{L}_{0}=\left(\partial_{\mu} \alpha\right) J_{\mu}+\left(\partial_{\mu} \alpha\right)^{2} \phi^{\star} \phi\tag{8.59}.$$

My question is, why the variation $\delta\mathcal L_0$ equals to $(\partial_\mu\alpha) J_\mu$ in the first order, and how you could obtain higher order corrections.
 A: This answer is going to focus on the concepts behind the equation you wrote down for the current, not the specific $U(1)$ example you give.
General variation of Lagrangian
Let's first consider a general transformation, parameterized by $\delta \lambda$ (I am intentionally going to save the symbol $\alpha$ for later). If we have $N$ fields $\phi_n$, where $n=1, 2, \cdots, N$, then under this transformation, the fields transform in the following way, to linear order in $\delta \lambda$
\begin{equation}
\delta \phi_n = \frac{\delta \phi_n}{\delta \lambda} \delta \lambda
\end{equation}
Meanwhile, the Lagrangian transforms as
\begin{eqnarray}
\delta \mathcal{L} &=& \sum_n \frac{\partial \mathcal{L}}{\partial \phi_n}\frac{\delta \phi_n}{\delta \lambda} \delta \lambda + \frac{\partial \mathcal{L}}{\partial \partial_\mu \phi_n} \partial_\mu \left( \frac{\delta \phi_n}{\delta \lambda} \delta \lambda \right) 
\end{eqnarray}
If we wanted to derive the Euler-Lagrange equations, at this stage we would integrate by parts and drop a boundary term on the second term. However, in this answer, we will not want to do that. Instead, note that we can use the product rule on the last term (this will come back later)
\begin{equation}
\partial_\mu \left( \frac{\delta \phi_n}{\delta \lambda} \delta \lambda \right) = \partial_\mu \left( \frac{\delta \phi_n}{\delta \lambda} \right) \delta \lambda +   \frac{\delta \phi_n}{\delta \lambda} \partial_\mu \delta \lambda  
\end{equation}
A clever derivation of the Noether current
Now, consider a continuous global symmetry parameterized by $\alpha$. A global symmetry means (a) that the parameter $\alpha$ does not depend on spacetime (so $\partial_\mu \alpha=0$) -- this is the "global" part, and (b) that $\delta \mathcal{L}=0$ (up to a total derivative$^\star$) under a transformation by $\alpha$ -- this is the "continuous symmetry" part. Note that we did not use the Euler-Lagrange equations to conclude that $\delta \mathcal{L}$ vanishes up to a total derivative. Instead, in this context, $\delta {\mathcal{L}}=0$ follows from the fact that $\alpha$ is associated with a symmetry.
Now, let us promote the parameter $\alpha$ of the global symmetry to a function $\alpha(x)$. The Lagrangian will generally not be invariant under this transformation. What can we say about the variation of $\mathcal{L}$ under $\alpha(x)$? Well, it must be the case that (up to a total derivative$^\star$) the variation is proportional to $\partial_\mu \alpha(x)$, because the variation must vanish when $\partial_\mu \alpha(x)=0$ as a consequence of the global symmetry. Furthermore, by Lorentz invariance, when $\partial_\mu \alpha$ appears in $\mathcal{\delta L}$, it must be contracted with a Lorentz vector. Let's call this vector $J^\mu$ (so far, this is not a current, but you can probably see where this is headed!). In equations, we can express this as
\begin{equation}
\delta \mathcal{L} = J^\mu \partial_\mu \alpha(x)
\end{equation}
Now we can integrate this by parts to get the variation
\begin{equation}
\delta \mathcal{L} = - \partial_\mu J^\mu \alpha(x)
\end{equation}
Now, because $\delta \mathcal{L}=0$ for any variation when the Euler-Lagrange equations are satisfied, and since $\alpha(x)$ is just a special case of a variation, it must be the case that $\partial_\mu J^\mu = 0$ when the Euler-Lagrange equations are satisfied. Thus, we identify $J^\mu$ as the conserved current in Noether's theorem.

$^\star$ While it's not important in your $U(1)$ example, so I won't pay much attention to it in this answer, this total derivative term is sometimes important to account for when computing the Noether current.
Computing $J^\mu$
As a final step, we can put things together and give the explicit formula for $J^\mu$, using our formula for a general variation $\delta \lambda$. Setting $\delta \lambda=\alpha(x)$, we have
\begin{eqnarray}
\delta \mathcal{L} &=& \sum_n \frac{\partial \mathcal{L}}{\partial \phi_n}\frac{\delta \phi_n}{\delta \alpha(x)}\alpha(x) + \frac{\partial \mathcal{L}}{\partial \partial_\mu \phi_n} \partial_\mu \left[\frac{\delta \phi_n}{\delta \alpha(x)} \alpha(x) \right] \\
&=& (...) \alpha(x) + \sum_n \frac{\partial \mathcal{L}}{\partial \partial_\mu \phi_n} \frac{\delta \phi_n}{\delta \alpha(x)} \partial_\mu \alpha(x)
\end{eqnarray}
the term in brackets $(...)$ multiplying $\alpha(x)$ must vanish identically (not using the Euler-Lagrange equations) as a consequence of the global symmetry, because of the general argument we gave above. (You can check that this cancellation happens explicitly in your massive complex scalar field example -- it is a direct consequence of the global $U(1)$ symmetry in that case).
Comparing this expression with the formula $\delta \mathcal{L} =  J^\mu \partial_\mu \alpha(x)$, we identify
\begin{equation}
J^\mu = \sum_n \frac{\partial \mathcal{L}}{\partial \partial_\mu \phi_n} \frac{\delta \phi_n}{\delta \alpha(x)} 
\end{equation}
as desired.
A: 
$$\mathcal{L}_{0}=\left(\partial_{\mu} \phi\right)^{\star}\left(\partial_{\mu} \phi\right)-m^{2} \phi^{\star} \phi$$


and


$$\delta \mathcal{L}_{0}=\left(\partial_{\mu} \alpha\right) J_{\mu}+\left(\partial_{\mu} \alpha\right)^{2} \phi^{\star} \phi\tag{8.59}.$$


I'm not entirely sure how these two equations follows, my best guess is that one can identify $\left(\partial_{\mu} \alpha\right) J_{\mu}$ as $\delta \alpha$...

No, $\left(\partial_{\mu} \alpha\right) J_{\mu}$ is not $\delta\alpha$; it says in the book and in your post that $\left(\partial_{\mu} \alpha\right) J_{\mu}$ is the first order transformation of the Lagrange density ($\delta\mathcal L_0$).
Regarding the first equation above, it also says in your book that $\mathcal L_0$ is the scalar QED Lagrange density for $A=0$ (or equivalently $e=0$). This is also something you wrote in your post.
Regarding the second equation you asked about, Schwartz is literally just making the following replacement in $\mathcal L_0$
$$
\phi(x) \to e^{-i\alpha(x)}\phi(x)
$$
and expanding the terms out...
Just use the fact that
$$
J_\mu = i\left( \phi\partial_\mu\phi^* - \phi^*\partial_\mu\phi \right)
$$
