I'm trying to show that for coordinates $x:U \to \mathbb{R}^n$ centered at $p \in (M,g)$, the property $g_{ij}(0)x^j = g_{ij}(x)x^j$ characterizes normal coordinates - i.e. straight lines through the origin are geodesics.
However, I don't know what $g_{ij}(0)x^j$ even means - is it just the matrix of coefficients of $g$ multiplied by some vector? Does it have any geometric meaning? $g$ is a map from $TM$ so it's weird that it's being applied to a point in the image of $M$ in $\mathbb{R}^n$. I get that $T_pM$ is diffeomorphic to some neighborhood of $p$, but I'm struggling to see the context for this property (it was shown in class and I wasn't in class).
If the proof is just symbol manipulation - something like taking the derivative of $g_{ij}(c(t))c^j(t)$ where $c(t)$ is a parametrized line in the direction of some vector $c$, then that's fine, I'm just wondering if there's more geometric context in the property/condition.
Here's one way to keep it straight. One can identify $T_p M$ with $\mathbb R^n$ through the correspondence $p \mapsto \mathbf 0$ and $v^i \frac{\partial}{\partial x^i} \mapsto v^i \mathbf e_i,$ where $\mathbf e_i$ as usual is the standard basis vector in the $i$th direction. In other words, $\mathbb R^n$ is a collection of vectors with which you may differentiate functions on $M$ at $p$. In this way, you can think of your coordinates $x: U \to \mathbb R^n$ as being from $U \to T_p M$.
From this viewpoint, the best way to think about normal coordinates - coordinates that map geodesics issuing from $p$ to straight lines issuing from $\mathbf 0$ - are coordinates such that the metric components have the special form $$ g_{ij} dx^idx^j = dr^2 + h_{ij} dS^{n-1},$$ where $h_{ij}$ is some metric on the unit sphere. A fancy way of saying this, if you know the lingo, is that the exponential map is a radial isometry. Why does this show that Euclidean straight lines get mapped to geodesics? The proof is in most introductory Riemannian Geometry texts, including Peter Petersen's Riemannian Geometry, but the basic idea is that you may verify that a straight line segment from $\mathbf 0$ to any point in $\mathbb R^n$ satisfies the geodesic equation w.r.t. the metric $g_{ij}$. On the other hand, if these straight euclidean lines are known to satisfy the geodesic equation, one may deduce the form of the metric above. In this way, normal coordinates are the analogue of polar coordinates in $\mathbb R^n$.
Now I can give you your characterization. First, $g_{ij}(0)$ is a symmetric, positive definite matrix so there exists a linear change of coordinates such that $g_{ij}(0) = \delta_{ij},$ the identity matrix. (Proof: From the spectral theorem of linear algebra we get $g_{ij}(0) = Q\Lambda Q^{-1},$ where $\Lambda$ is a diagonal matrix with positive entries $\lambda_i$. Define $\Lambda^{-1/2}$ to be the matrix with diagonal entries $\lambda_i^{-1/2}$. Then the linear coordinate change is $ y = Q\Lambda^{-1/2}x$.) Therefore, without loss of generality, modify $x$ by this coordinate change and assume $g_{ij}(0) = \delta_{ij}$.
Suppose first that $g_{ij}(x)x^j = g_{ij}(0)x^j = \delta_{ij}x^j = x_i.$ I will first note that for any point $x \in \mathbb R^n$, we may write $x = (x^1, \ldots, x^n) = \sum_i x^i \mathbf e_i = x^i.$ This is the meaning of $x^i.$ We may also think of $x$ as a vector from $\mathbf 0$ to the point $x$. This is, of course, also the radial vector $\mathbf r = x^i$. I will write $r = || x||$ for the length of $\mathbf r$. Now $g_{ij}(x)x^j$ is the function $g_x(\mathbf r, \cdot) = g_x(x^i, \cdot)$, which by our assumption is the function $g_p(\mathbf r, \cdot) = \delta_{ij}x^j = x_i$. If you are familiar with the more formal language, $x_i$ is called a co-vector. Write $\hat{\mathbf r}$ for the unit vector $\frac{\mathbf r}{r}$. Then for any $x \in \mathbb R^n$, we have $$g_x(\hat{\mathbf r},\hat{\mathbf r}) = \frac{1}{r^2} g_x(\mathbf r, \mathbf r) = \frac{1}{r^2}x_i\hat{\mathbf r} = \frac{1}{r^2}x_ix^i = \frac{r^2}{r^2} = 1,$$
which shows that $$g_{ij}(x) = dr^2 + \textrm{ other terms. }$$ Now suppose that $\mathbf w= w^i \mathbf e_i = w^i$ is any vector that is Euclidean-orthogonal to $\mathbf r,$ so that $\mathbf r \cdot \mathbf w = \delta_{ij} x^j w^i = 0$. Then we may calculate that $$g_x(\mathbf r, \mathbf w) = \delta_{ij}x^j w^i = 0,$$ which shows that Euclidean-orthogonal vectors to $\mathbf r$ are also $g_{ij}$-orthogonal to $\mathbf r$. This says that $$ g_{ij}(x) = dr^2 + h_{ij}dS^{n-1},$$ for some metric $h_{ij}$ on the unit sphere $S^{n-1}$. In other words, we have shown that there are no `interaction' terms of the form $drdw$.
Your intuition, therefore, is that the expression $g_{ij}(x)x^j = g_{ij}(0)x^j$ is telling you that the metric stays constant in the radial direction. Since every $g_{ij}(0)$ is the Euclidean metric after a linear change of coordinates, this means that the expression $g_{ij}$ `stays Euclidean' in the radial direction: $g$ is a radial isometry. The intuition is clear and beautiful.
I could also go the other direction and show you that $g$ in normal coordinates $g_{ij} = dr^2 + h_{ij}S^{n-1}$ implies that $g_{ij}(x)x^j = g_{ij}(0)x^j.$ However: It is late in the evening and I am getting tired. I would be happy to update with more info later if you are interested.