Why is the second derivative of $f:\mathbb R^n \to \mathbb R$ a function $\mathbb R^n \times \mathbb R^n \to \mathbb R$? In my multivariable analysis class notes my teacher wrote the following:

Let $D \subseteq \Bbb R^n$, $p \in D$ and $f:D \to \mathbb R$ be a class $C^2$ function. The function: $$\mathbb R^n \times \mathbb R^n \to \mathbb R$$ $$((y_1,...,y_n),(z_1,...,z_n)) \mapsto \sum_{i,j = 1}^n y_i z_j \frac{\partial ^2 f}{\partial x_j\partial x_i}(p)$$
Is the second derivative of the function $f$ at the point $p$, denoted ad $f''(p)$.

Why is the second derivative a function with domain $\mathbb R^n \times \mathbb R^n$?
If $\mathcal L(A,B)$ denotes the set of all bounded linear operators from $A$ to $B$, shouldn't the second derivative of $f$ at $p$ be a function $f'' : \mathbb R^n \to \mathcal{L}( \mathbb R^n,\mathcal{L}(\mathbb R^n, \mathbb R) )$ instead?
 A: Firstly, you are right in saying that the second differential $f''(p)$ must be, by definition, a linear map $f''(p)\colon \mathbb{R}^{n}\to \mathcal{L}(\mathbb{R}^{n},\mathcal{L}(\mathbb{R}^{n},\mathbb{R}))$. Nevertheless, we have a way to identify the target space with the set of all bilinear forms over $\mathbb{R}^{n}$ (which would seem to be simpler to understand).
We need a bit of abstract nonsense from Linear Algebra. Suppose $V$ is a finite dimensional space over, say, the real numbers. A bilinear form $g\colon V\times V\to \mathbb{R}$ is a function which is essentially linear in both variables. In other words, fixed $v\in V$, the functions $g(\cdot,v)$ and $g(v,\cdot)$ must be linear forms.
Suppose $g$ is a bilinear form over $V$. Because $g(v,\cdot)$ is linear for every $v\in V$, the map $v\mapsto g(v,\cdot)$ is a well defined map from $V$ to $\mathcal{L}(V,\mathbb{R})$. Further, because $g(\cdot,w)$ is linear for any $w$, this means that the map $v\in V\mapsto g(v,\cdot)\in \mathcal{L}(V,\mathbb{R})$ is a linear map.
In essence, we have constructed a map $\Phi\colon\mathcal{L}^{2}(V,\mathbb{R})\to \mathcal{L}(V,\mathcal{L}(V,\mathbb{R}))$ given as follows:
$$
g\in\mathcal{L}^{2}(V,\mathbb{R})\mapsto\Phi(g)\in \mathcal{L}(V,\mathcal{L}(V,\mathbb{R})), \quad \Phi(g)(v)=g(v,\cdot).
$$
Try to check that $\Phi$ is a linear map (the notation might seem a little overpowered, but it's all a matter of unwinding the definitions). The key here is that $\Phi$ is, in reality, an isomorphism of vector spaces. Indeed, the map $\Lambda\colon \mathcal{L}(V,\mathcal{L}(V,\mathbb{R}))\to \mathcal{L}^{2}(V,\mathbb{R})$ given by
$$
h\in \mathcal{L}(V,\mathcal{L}(V,\mathbb{R}))\mapsto \Lambda(h) \in \mathcal{L}^{2}(V,\mathbb{R}), \quad \Lambda(h)(v,w)=h(v)(w)
$$
is the inverse of $\Phi$ (try to check that $\Lambda$ is well defined and the inverse of $\Phi$, you don't need to prove linearity).
So, since $\Phi$ is a linear isomorphism, we can identify $\mathcal{L}^{2}(V,\mathbb{R})$ and $\mathcal{L}(V,\mathcal{L}(V,\mathbb{R}))$. We usually say that $\Phi$ is a natural isomorphism, because to construct it we didn't need more information about $V$ except that it is a vector space (in other words, what we just did works for any generic vector space over any field).
Let's return to our original question. Suppose we have a $\mathscr{C}^{2}$ map $f\colon D\subseteq \mathbb{R}^{n}\to \mathbb{R}$, and let $p\in D$ (I'm assuming $D$ is open). Since the first derivative is a differentiable function $f'\colon D\to \mathcal{L}(\mathbb{R}^{n},\mathbb{R})$, the second derivative $f''(p)$ lives in $\mathcal{L}(V,\mathcal{L}(V,\mathbb{R}))$. Nevertheless, because of the isomorphism we wrote down earlier, we can identify $\mathcal{L}(V,\mathcal{L}(V,\mathbb{R}))$ with $\mathcal{L}^{2}(V,\mathbb{R})$, so knowing $f''(p)$ is the same as knowing $\Phi(f''(p))$.
Now, let's think about what the bilinear form $\Phi(f''(p))$ looks like. From Linear Algebra, we know that it must admit an expression of the form
$$
\Phi(f''(p))((x_{1},...,x_{n}),(y_{1},...,y_{n}))=\sum_{i,j=1}^{n}a_{ij}x_{i}y_{j}
$$
where $a_{ij}=\Phi(f''(p))(e_{i},e_{j})=f''(p)(e_{i})(e_{j})$ (here, $\{e_{i}\}$ is the standard basis). Let us compute $a_{ij}$ by unwinding the definitions a little.
We know, since the differential at a vector is equal to the directional derivative, that
$$
f''(p)(e_{i})(e_{j})=\dfrac{d}{dt}\Big|_{t=0}f'(p+te_{i})(e_{j})=\dfrac{d}{dt}\Big|_{t=0}\dfrac{\partial f}{\partial x_{j}}(p+te_{i})=\dfrac{\partial^{2}f}{\partial x_{i}\partial x_{j}}(p).
$$
This means that the second differential $f''(p)$ is identified with the bilinear form
$$
\Phi(f''(p))\colon (x,y)\in \mathbb{R}^{n}\times \mathbb{R}^{n}\mapsto \sum_{i,j=1}^{n}\dfrac{\partial^{2}f}{\partial x_{i} \partial x_{j}}(p)x_{i}y_{j}.
$$
Once you are aware of the identification given by $\Phi$, you just omit it and say that $f''(p)$ is the above map (and you're not losing information since $\Phi$ is an isomorphism!). Similarly, for higher order differentials, you would identify $f^{(k)}(p)$ with a $k$-multilinear form, which is far easier than working with iterated linear map spaces.
Hope this helps!
