We denote by $e_f(x) \geq 1$ the order of vanishing of a cusp form $f$ at a cusp $x$.
There is a newform $f$ of weight 2 and level $\Gamma_0(625)$ with coefficients in $\mathbb{Q}(\sqrt{5})$ and Fourier expansion
\begin{equation*}
f = q + \left(\frac{-1-\sqrt{5}}{2}\right) q^2 + \left(\frac{-3-\sqrt{5}}{2}\right) q^3 + \left(\frac{-1+\sqrt{5}}{2}\right) q^4 \ldots
\end{equation*}
One can show that
\begin{equation*}
e_f(a/25) = \begin{cases} >1 & \textrm{if } a \equiv \pm 1 \textrm{ mod } 5 \\
1 & \textrm{if } a \equiv \pm 2 \textrm{ mod } 5.
\end{cases}
\end{equation*}
The cusps $a/25$ are Galois conjugate. If we denote by $\sigma$ the non-trivial automorphism of $\mathbb{Q}(\sqrt{5})$, we get $e_{f^\sigma}(1/25)=1$ while $e_f(1/25)>1$. Note here that $\mathbb{Q}(\sqrt{5})$ is contained in the cyclotomic field $\mathbb{Q}(\zeta_{625})$.
I should add that determining the order of vanishing of a modular form at an arbitrary cusp is a difficult question in general. You can certainly guess it numerically by just estimating the decay rate of the modular form at the given cusp. However, proving it rigorously is much more delicate. For the modular form above this involves computations with the local automorphic representation associated to $f$, which in this case is supercuspidal and is described by an irreducible representation of $\mathrm{GL}_2(\mathbb{Z}/25\mathbb{Z})$ of dimension 20 with coefficients in $\mathbb{Q}(\sqrt{5})$.
For other examples, one can show that if $f$ is a newform on $\Gamma_0(p^4)$ with $p \geq 5$ prime, such that the local automorphic representation $\pi_{f,p}$ is a (ramified) principal series, then $e_f(a/p^2)$ is equal to 1 for approximately half of the $a$ in $(\mathbb{Z}/p^2\mathbb{Z})^\times$, while it is $>1$ for the other half. This follows from the vanishing of certain character sums modulo $p^2$ proved by Paul Nelson (unpublished). So we get further examples by looking at newforms of conductor divisible by a high power of a prime. The examples of level 567 and 891 in my comment are of this kind.
There are algorithms to compute the Fourier expansion of a modular form at arbitrary cusps (see this MO question), but they are either numerical or may become slow when the conductor gets large.
For theoretical results, you can look at the article of Corbett and Saha, On the order of vanishing of newforms at cusps.