After playing around with power series $S_{n, p}=\sum_{k=0}^nk^p$ and finding a formula for each natural power, I tried to also find formulae for imaginary power but couldn't figure out how.
$S_n=\sum_{k=0}^nk^i$, where $i=\sqrt{-1}$
I tried expressing the series using cis and got $S_n=\sum_{k=0}^ncis(ln(n))$ but this didn't help me at all. Any hints or full answer will be appreciated.
By the Euler-Maclaurin formula,
$$\sum_{k=1}^nk^i=\zeta(-i)+\frac{n^{i+1}}{i+1}+\frac{n^i}2+\frac{in^{i-1}}{12}-\frac{i(i-1)(i-2)n^{i-3}}{720}+\mathcal O(n^{i-5})$$
this is not an exact formula, but rather an approximation to the problem, where $\zeta(-i)\approx0.0033002237+0.4181554491i$ is the Riemann zeta function. A few values are given to show how it approximates:
$$\begin{array}{c|c}n&\sum&\zeta\\\hline1&1&0.49913+0.00010i\\2&1.76924+0.63896i&1.38459+0.31946i\\3&2.22407+1.52954i&1.99665+1.08425i\end{array}$$
Ok, not so good for small values, but it gets better.
$$\begin{array}{c|c}n&\sum&\zeta\\\hline10&0.41898+7.84548i&0.37600+7.47349i\\20&-8.93237+11.8340i&-8.43768+11.76128i\\30&-18.82708+10.93402i&-18.34384+11.06237i\end{array}$$
As you can see, the accuracy increases as $n\to\infty$.
Alternative forms include
$$\sum_{k=1}^nk^i=\zeta(-i)-\zeta(-i,n+1)=H_n^{(-i)}$$
where $\zeta(s,q)$ is the Hurwitz zeta function and $H_n^{(p)}$ is the generalized harmonic number.
It has been pointed out that OP's question is about the finite sum, so my answer does not address the question.
The sum $\sum^\infty_{k=1} k^{-s}$ is absolutely convergent only when $\Re s>1$. So the infinite sum $\sum^\infty k^i$ does not converge. However in the region $\Re s>1$, the sum is analytic, so admits a unique continuation to a maximal domain on the complex domain. This continuation is called the Riemann zeta function, and its value at -i is in some sense related to the infinite sum $\sum^\infty k^i$. That value is approximately 0.0033 + 0.418$i$.
By the way, when $s$ is a natural number, the formula for the sum is called Faulhaber's formula.
$$ \sum_{k=1}^nk^s=\frac{1}{s+1}\sum_{j=0}^s(-1)^j{s+1\choose j}B_jn^{s+1-j} $$ where $B_j$ is the $j$th Bernoulli number. This formula does not make any sense for values of $s$ which are not natural numbers, though.
