How to prove the raising/lowering indices operation?
Also I've searched for it in books like Carroll's or Lawden's, but it's given pretty much as if it would be a definition.
Because it is. No need for differential geometry, linear algebra is sufficient here:
At a given point of space-time, the tangent space is just a vector space, the cotangent space its dual (ie the space of real-valued linear functions).
Given a basis $\{{\bf b}_i\}$ of the vector space, we define the dual basis $\{\beta^i\}$ via $$ \beta^i({\bf b}_j)=\delta^i_j $$ and can expand any vector and dual vector in terms of these (and implicitly do so in Einstein notation, using the basis induced by a coordinate chart).
As the vector spaces have the same (finite) dimension, they are of course isomorphic, eg via the isomorphism $$ {\bf b}_i\mapsto\beta^i $$ The problem is that this is not a canonical isomorphism, ie it depends on the choice of basis: Use a different one, and you'll generally end up with a different dual vector.
So we need another way to specify an isomorphism. Using a non-degenerate bilinear form like the metric tensor to 'lower' indices is such a mapping from vectors to dual vectors, which by definition comes with an inverse ('raising' indices).
As suggested in the comments, "lowering an index" is just coordinate notation for the isomorphism $\flat:TM\to T^*M$ between defined by $$X^\flat (Y) = \langle X,Y\rangle\text{,}$$ where $X$ and $Y$ are arbitrary vectors.
I've tried using the definition of the metric ${g_{\alpha\beta}=\hat{\mathbf{e}}_\alpha\cdot\hat{\mathbf{e}}_\beta}$ where $\{\hat{\mathbf{e}}_\mu\}$ is a basis of the space, ...
That's all you need to show that the above is just the same as index-lowering. You can expand the vectors in an arbitrary basis and use the fact that the metric tensor components are precisely the inner products of the coordinate vectors, as you say here.
Similarly, "raising an index" is just coordinate notation for the isomorphism $\sharp:T^*M\to TM$ defined by $$\omega(Y) = \langle\omega^\sharp,Y\rangle\text{,}$$ where $\omega$ is an arbitrary covector.
These are sometimes called musical isomorphisms, although this label is not particularly common in physics.
Notation
$C(M)$ represents smooth functions $f:M\rightarrow \mathbb{R}$. $\tau$ represents a tensor field, and $\mathcal{T}^r_s(M)$ represents the set of all $(r,s)$ tensor fields on $M$. The function $\langle \cdot,\cdot \rangle$ represents the operation of the metric $g(\cdot,\cdot)$. And denote $g=\sum_{i,j} g_{ij} dx^i\otimes dx^j$ where $\partial_i$ and $dx^i$ are the respective coordinate vector fields and coordinate one-forms. The symbols $\mathfrak{X}(M)$ and $\mathfrak{X}^*(M)$ represent the modules of smooth vector fields and one-forms on $M$ respectively.
Proposition
Let $(M,g)$ be a pseudo-Riemannian manifold. If $X\in\mathfrak{X}(M)$, then let $\theta_X\in \mathfrak{X}^*(M)$ such that $\theta_X(Y)=\langle X,Y \rangle \ \ \forall Y\in \mathfrak{X}(M)$. The function $\phi:\mathfrak{X}(M)\rightarrow \mathfrak{X}^{*}(M)$ mapping $X\mapsto \theta_X$ is a $C(M)$-linear isomorphism. Given such an isomorphism, we say that $X$ and $\theta_X$ are metrically equivalent.
Theory
This implies that one can change a vector field to a unique one-form on a pseudo-Riemannian manifold since the corresponding pairs $X\leftrightarrow \theta_X$ contain the same information up to isomorphism, which implies $\mathcal{T}_0^1(M)\simeq \mathcal{T}_1^0(M)$. Consider approaching this idea to arbitrary dimensions of tensor fields. Let $a\in \{1,\ldots,r\}$ and $b\in\{1,\ldots,s\}$, let $\tau\in \mathcal{T}^r_s(M)$, then express the corresponding lowered index tensor field $^a_{b*}\tau\in \mathcal{T}^{r-1}_{s+1}(M)$ as \begin{align*} (^{a}_{b*} \tau ) (\theta ^{1} ,\dotsc ,\theta ^{r-1} ,X_{1} ,\dotsc ,X_{s+1} )\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \ \ \ \ \ \ \ \ \ \ \ \ \ =\tau (\theta ^{1} ,\dotsc ,\theta ^{a-1} ,X^{*}_{b} ,\theta ^{a+1} \dotsc ,\theta ^{r-1} ,X_{1} ,\dotsc ,X_{b+1} ,X_{b-1} ,\dotsc ,X_{s} ) \end{align*} where $X_b\leftrightarrow X_b^*$ is the metrically equivalent pairs. The function $\psi(\tau)=~^a_{b*}\tau$ essentially transform a vector field in the $b^{\text{th}}$ position to a one-form, and placing it in the $a^{\text{th}}$ position. The function is clearly well defined since there $r$ number of one-form arguments and $s$ number of vector field arguments on the right hand side. Given that $\psi$ is an isomorphism, the inverse $\psi^{-1}$ clearly exists, and it represents the raising index operation. To be precise, denote the raising index function as $\psi^{-1}(\tau)= ~ ^{a*}_b\tau$ where $^{a*}_b\tau\in \mathcal{T}^{r+1}_{s-1}$.
Derivation
However, we have not yet describe the computation which coverts a vector field $X_b$ to its metrically equivalent one-form $X_b^*$. With respect to a coordinate chart $(U,(x^i))$, it is given that $dx^i(\partial_j)=\partial_j(x^i)=\delta^i_j$. For any $Y\in \mathfrak{X}(M)$, we have $$\langle \partial_k,Y \rangle = \sum_{i,j} g_{ij} dx^i(\partial_k)dx^j(Y) = \sum_j g_{kj} dx^j (Y) \ \ \forall k\in \{1,\ldots,n\}$$ by using the preceding proposition, the unique isomorphism can be defined as $$\psi:\partial_k\mapsto \sum_j g_{kj}dx^j$$ with the corresponding inverse $$\psi^{-1}:dx^k\mapsto \sum_j g^{kj}\partial_j$$ Therefore, any vector field $X\in \mathfrak{X}(M)$ can be "converted" by $$\psi(X)=\psi(\sum_i X^i \partial_i)=\sum_i X^i\psi(\partial_i)= \sum_{i,j} X^i g_{ij}dx^j$$ which uniquely defines a one-form in $\mathfrak{X}^*(M)$. The inverse can be computed for a one-form $\theta\in\mathfrak{X}^*(M)$ similarly to obtain $$\psi^{-1}(\theta)=\sum_{i,j}\theta_ig^{ij}\partial_j$$