Sign of the totally anti-symmetric Levi-Civita tensor $\varepsilon^{\mu_1 \ldots}$ when raising indices

In flat spacetime, the isomorphism between contra- and covariant components is furnished by the Minkowski metric $\eta_{\mu\nu}$. The Minkowski metric in $n$ spacetime dimensions is just $\operatorname{diag}(-1,1,\dots,1)$. Thus for all dimensions, $$\operatorname{det}\eta=-1$$ For any $n\times n$ matrix $A_{\mu\nu}$, there is a well-known theorem $$\epsilon_{\mu_1\cdots\mu_n}\operatorname{det}A=\sum_{i}\sum_{\nu_i}\epsilon_{\nu_1\cdots\nu_n}A_{\nu_1\mu_1}\cdots A_{\nu_n\mu_n}$$ where we have made no distinction between covariance and contravariance. But suppose $A=\eta$. Then all the multiplications by $\eta$ will raise the indices on $\epsilon$. And since $\operatorname{det}\eta=-1$ in any spacetime dimension, we get $$\epsilon^{\cdots}=-\epsilon_{\cdots}$$ where the dots represent any number of indices.

Sean Carroll's Spacetime and Geometry has a thorough discussion of this, and, even better, this discussion is in the lecture notes that turned into the book (see Chapter 2: Manifolds).

In full generality (at least for any right-handed coordinate system), we start off with the symbol $\tilde{\epsilon}_{\mu_1\cdots\mu_n} \equiv [\mu_1\,\cdots\,\mu_n]$, which is $0$ or $\pm1$ depending on the sign of $\mu_1\,\cdots\,\mu_n$ as a permutation of $0\,\cdots\,(n-1)$. Then the tensor with lower indices obeys $$ \epsilon_{\mu_1\cdots\mu_n} = \sqrt{\lvert g \rvert} \tilde{\epsilon}_{\mu_1\cdots\mu_n} = \sqrt{\lvert g \rvert}\ [\mu_1\,\cdots\,\mu_n], $$ where $g$ is the determinant of the metric.

Some people define a symbol with upper indices that is the same as that with lower indices, but on the other hand some people include an extra factor of $\operatorname{sgn}(g)$. In fact, this is ambiguous enough that Carroll states one version in the linked notes, and the other version in the final published book.

In any event, the tensor with upper indices is unambiguously $$ \epsilon^{\mu_1\cdots\mu_n} = \frac{1}{g} \epsilon_{\mu_1\cdots\mu_n} = \frac{\operatorname{sgn}(g)}{\sqrt{\lvert g \rvert}} [\mu_1\,\cdots\,\mu_n]. $$

In the end, you should be using the full tensors in calculations, at least if everything else in the formula is a true tensor. In flat spacetime, the metric is $\operatorname{diag}(-1,+1,+1,+1)$, so $g = -1$ and we have \begin{align} \epsilon_{\mu_1\cdots\mu_n} & = \phantom{-}[\mu_1\,\cdots\,\mu_n], \\ \epsilon^{\mu_1\cdots\mu_n} & = -[\mu_1\,\cdots\,\mu_n]. \end{align}