The set $V= \{ A\in M_n(\mathbb{C})\ |\ A\bar A = I\}$ is a smooth submanifold of $M_n(\mathbb{C})$ with real dimension $n^2$.

Proof: Consider the involution $\iota:\mathrm{GL}(n,\mathbb{C})\to \mathrm{GL}(n,\mathbb{C})$ defined by
$$
\iota(A) = (\bar A)^{-1}.
$$
This is an anti-holomorphic involution of the complex manifold $\mathrm{GL}(n,\mathbb{C})$, and its fixed locus is precisely $V$. Thus, $V$ is a totally real submanifold of $\mathrm{GL}(n,\mathbb{C})$ with (real) dimension $n^2$.

Also: While it's true that the map $f:M_n(\mathbb{R})\to \mathrm{GL}(n,\mathbb{C})$ defined by
$$
f(a) = \exp(ia)
$$
has its image in $V$, it is not a 'parametrization' everywhere, i.e., $f$ is not a local diffeomorphism everywhere. While this is true on a neighborhood of $0\in M_n(\mathbb{R})$, at other places, the map $f$ definitely is not a local diffeomorphism. For example, if $a$ has $n$ eigenvalues of the form $2k_i\pi$ (where $k_1,\ldots, k_n$ are integers, not all zero), then $f(a) = I$, but the space of such real matrices $a$ has positive dimension for any $n$-tuple $(k_1,\ldots, k_n)$ for which not all of the $k_i$ are equal.

It turns out that $f$ is surjective. The proof uses the fact that, for $A\in V$, we have $A\bar A = \bar A A = I$, so, in particular, $A$ and $\bar A$ commute and hence can be put simultaneously in Jordan normal form by a real conjugation. Then, breaking $\mathbb{C}^n$ into a sum of complexifications of real subspaces according to the eigenvalues of $A$ and using the semi-simple/nilpotent decomposition appropriately, one can reduce to dealing with the upper triangular case, and the result can then be proved using simple facts about power series in commuting nilpotent variables. Details upon request (**see below**).

Of course, all of this is probably a special case of known facts about affine symmetric spaces. I have now realized that $V$ is simply the Cartan embedding for the affine symmetric space $\mathrm{GL}(n,\mathbb{C})/\mathrm{GL}(n,\mathbb{R})$. This Cartan embedding
$$
\sigma:\mathrm{GL}(n,\mathbb{C})/\mathrm{GL}(n,\mathbb{R})\longrightarrow V\subset
\mathrm{GL}(n,\mathbb{C})
$$
is given by $\sigma(B\cdot\mathrm{GL}(n,\mathbb{R})) = B\,(\,\overline B\,)^{-1}$. The exponential map we have been discussing is just the geodesic mapping of this affine symmetric space, so, probably all of this follows from general theory.

**Details:** Let $A\in\mathrm{GL}(n,\mathbb{C})$ satisfy $A\overline{A} = I$. Let $\lambda\in\mathbb{C}$ be an eigenvalue of $A$ of multiplicity $m\le n$ and let
$$
V_\lambda = \{ v\in\mathbb{C}^n \ |\ (A{-}\lambda)^mv = 0\ \}.
$$
be the associated generalized eigenspace. Of course, $\mathbb{C}^n$ is the direct sum of these generalized eigenspaces. Since $\overline{A}$ commutes with $A$, it preserves each of these subspaces. Moreover, one clearly has
$$
\overline{V_\lambda} = V_{1/\bar\lambda},
$$
so each of the spaces $V_\lambda+V_{1/\bar\lambda}$ is invariant under conjugation and hence is the complexification of a real subspace $W_\lambda\subset \mathbb{R}^n$. It follows that, by conjugating by a real matrix, we can assume that $A$ (and hence $\overline A$) is in block diagonal form, so it suffices to treat the case where either $A$ has a single eigenvalue $\lambda$ satisfying $\lambda\bar\lambda=1$, so that $\mathbb{C}^n=V_\lambda$, or else $A$ has an eigenvalue $\lambda$ satisfying $\lambda\bar\lambda > 1$ and $\mathbb{C}^n=V_\lambda\oplus V_{1/\bar\lambda}$.

In either case, we can assume that $\lambda = r\ge 1$ is real, since, if $\lambda = r e^{i\theta}$, we can replace $A$ by $e^{-i\theta}A$ and show that $e^{-i\theta}A$ is in the image of $f$ (since $I$ commutes with everything and $e^{i\theta} I = \exp(i\theta I)$.)

First, consider the case $\lambda = 1$. Then $A = C + i S$ where $C$ and $S$ are real matrices and $((C-I) + iS)^n = (A-I)^n = 0$. Moreover, $N=C-I$ and $S$ commute since $I = A\bar A = C^2+S^2 +i(SC-CS) = I + 2N + N^2 + S^2$. Note that $N$ and $S$ are nilpotent commuting matrices. Since they satisfy
$$
2N + N^2 = -S^2,
$$
it follows that $N = p(-S^2)$ where $p(t) = -\tfrac12 t^2 + \cdots $ is the (unique) power series (with real coefficients) that satisfies $2p(t) + p(t)^2 = -t^2$. Since $S$ is nilpotent, it follows that $N = p(-S^2)$ expresses $N$ canonically as a polynomial in $S$. Now let $q(t)$ be the (unique) power series with real coefficients that satisfies $\sin q(t) = t$. Note that we must have $p(-q(t)^2) = \cos(t)-1$. Now, because $S$ is nilpotent, we have $S = \sin q(S)$, where $q(S)$ is a real polynomial in $S$. Putting all of this together, we have
$$
A = I + p(-S^2) + iS = \cos(q(S)) + i\sin(q(S)) = \exp(iq(S)).
$$
Finally, before leaving this special case, let us note that, because $A$ satisfies its characteristic polynomial $(A-I)^n = 0$, it follows that $\overline A = A^{-1}$ is expressed as a universal polynomial in $A$ with real coefficients, so $iS = \tfrac12(A-\overline{A})$ is also expressed as a universal polynomial with real coefficients. Since $q(t)$ is an odd power series, it satisfies $iq(it) = f(t)$ where $f$ has real coefficients, so there is actually a formula of the form $A = \exp( i g_n(A))$ for all $A \in \mathrm{GL}(n,\mathbb{C})$ satisfying $A\bar A = I$ and $(A-I)^n=0$ for some universal polynomial $g_n(t)$ with real coefficients that also has the property that $g_n(A)$ is a real $n$-by-$n$ matrix for all such $A$. (This remark will be used below.)

Now, finally, let us assume that $A$ has eigenvalues $\lambda = e^t$ and $1/\lambda = e^{-t}$ for some $t>0$. Then $V_\lambda$ and $\overline{V_\lambda} = V_{1/\lambda}$ are disjoint, complementary complex subspaces of $\mathbb{C}^n$ and so we must have that there exists a matrix $Q\in M_n(\mathbb{R})$ such that
$$
V_\lambda = \{ v + iQv \ | \ v\in \mathbb{R}^n \}.
$$
Because $V_\lambda$ is closed under multiplication by $i$, it follows that $Q^2 = -I$. Setting
$$
S = \cosh t + i\sinh t\,Q = \exp(i t Q),
$$
we see that $Sw= e^t w$ for all $w \in V_\lambda$ and $Sw = e^{-t}w$ for all $w \in V_{1/\lambda}$, so $S$ is the semi-simple part of $A$. Hence $S$ can be written as $S = s_\lambda(A)$, where $s_\lambda(t)$ a polynomial in $t$ with real coefficients (that depend on $\lambda$). Thus, $\bar S = s_\lambda(\bar A)$ can be written as a universal polynomial in $A$, implying that $Q$ itself can be written as a universal polynomial in $A$ and hence, in particular, it commutes with $A$ and $\bar A$. Now, writing
$$
A = \exp(i t Q) B,
$$
we see that $B\bar B = I$ and that $B$ can be written as a polynomial in $A$. Moreover, the eigenvalues of $B$ are now all equal to $1$, so $(B-I)^n=0$, and so $B = \exp(i g_n(B))$ (as per above), where $g_n(B)$ is a polynomial in $A$, which, therefore, commutes with $Q$ (which is a polynomial in $A$). Finally, we have
$$
A = \exp(i t Q)\exp(i g_n(B)) = \exp(i (t Q+g_n(B))),
$$
as desired.