If $$\delta S = \int \sqrt g F[\phi] \delta \phi\tag{1}$$
Then is it natural to define the functional derivative as follows,
$$\frac{\delta S}{\delta \phi} = F[\phi].\tag{2}$$
In particular does this definition satisfy the commutativity of the functional derivatives.
I understand that this is not the standard definition of a functional derivative, but if I define this way this makes a certain calculation I am doing much easier to control. So to Clarify what I want to know is that if I define the 'functional derivative' this way then if
$$S = \int \sqrt{g} L[\phi, g_{\mu \nu}]\tag{3}$$ then
is the following true?
$$\frac{\delta}{\delta\phi} \frac{\delta}{\delta g_{\mu \nu}} S = \frac{\delta}{\delta g_{\mu \nu}} \frac{\delta}{\delta\phi} S \tag{4} $$
This is my procedure for computing the second functional derivative suppose
$$\frac{\delta}{\delta \phi}S = E[\phi, g_{\mu \nu}]\tag{5}$$
and
$$\frac{\delta}{\delta g_{\mu \nu}}S = E^{\mu \nu}[\phi, g_{\mu \nu}]\tag{6}$$
(equations of motion) Then to compute second functional derivative we write, eg.
$$E[\phi, g_{\mu \nu}](x) = \int \sqrt{g} d^4 y E(y) \hat \delta(x-y)\tag{7}$$
where $$\hat \delta (x-y) = \frac{\delta(x-y)}{\sqrt{g}}\tag{8} $$ is the generalised delta function. After this we can use the same definition to compute
$$\frac{\delta} {\delta g_{\mu \nu}} \int \sqrt{g} d^4 y E(y) \hat \delta(x-y)\tag{9} .$$ Of course now the result would involve delta functions.
Best Answer
Comments to the question (v4):
First a disclaimer. Note that even for a smooth local functional, the existence of a functional/variational derivatives is not guaranteed, but depends on appropriately chosen boundary conditions.
Calculus of variations, functional/variational derivatives, Frechet derivatives, Gateaux derivatives, etc, is a huge mathematical topic. If OP is looking for rigor, then we suggest to post the question on Math.SE or Mathoverflow.SE. In this answer we will just make some heuristic comments and ignore boundary terms.
Let there be given an $n$-dimensional pseudo-Riemannian manifold $(M,g)$.
We can consider a Dirac delta distribution $$\tag{A} \delta_M(x,y)~=~\frac{\delta^n(x-y)}{\sqrt{|g(x)|}}, \qquad g(x)~:=~\det(g_{\mu\nu}(x)),$$ on the manifold, cf. OP's eq. (8). Here $x$ and $y$ denote two points in a coordinate neighborhood $U\subseteq M$, [and to simplify notation, they also denote the coordinates of the two points].
The Dirac delta distribution (A) is independent of the choice of coordinates, i.e. transforms as a scalar.
Next we would like to define the functional/variational derivative wrt. fields $\phi^{\alpha}(x)$ [which could include the metric field $g_{\mu\nu}(x)$].
The functional derivative $\frac{\delta_M S}{\delta\phi^{\alpha}(x)}$ should satisfy $$\tag{B} \delta S~=~ \int \! d\mu(x)~\frac{\delta_M S}{\delta\phi^{\alpha}(x)}\delta\phi^{\alpha}(x), \qquad d\mu(x)~:=~\sqrt{|g(x)|}d^nx,$$ for appropriate infinitesimal variations $\delta\phi^{\alpha}(x)$, cf. OP's eq. (1).
In particular, we have the rule $$\tag{C} \frac{\delta_M \phi^{\beta}(y)}{\delta\phi^{\alpha}(x)} ~=~\delta^{\beta}_{\alpha}~\delta_M(x,y). $$ When $\phi^{\alpha}$ is the metric field $g_{\mu\nu}$, note this subtlety.
When $\phi^{\alpha}$ and/or $\phi^{\beta}$ are the metric field $g_{\mu\nu}$, then commutativity of functional derivatives are generally not expected, since the manifold notion (B) of functional differentiation depends on the metric, cf OP's eq. (4). In fact, it is easy to find counterexamples.
For the rest of this answer we will assume that $\phi^{\alpha}$ and $\phi^{\beta}$ are not the metric field $g_{\mu\nu}$. It is then expected that the functional derivatives commute $$\tag{D} [\frac{\delta_M}{\delta\phi^{\alpha}}, \frac{\delta_M}{\delta\phi^{\beta}}]~=~0. $$ See also e.g. Ref. 1.
Example. If the functional is of the form $$\tag{E}S~=~\int \! d\mu(x)~L(\phi(x),\partial\phi(x),x),$$ and we define $$\tag{F}\nabla^{(x)}_{\mu} ~:=~ \frac{1}{\sqrt{|g(x)|}}d^{(x)}_{\mu}\sqrt{|g(x)|} ,\qquad d^{(x)}_{\mu}~:=~\frac{d}{dx^{\mu}},$$ then the 1st functional derivative reads $$\frac{\delta_M S}{\delta\phi^{\alpha}(x)} ~=~\frac{\partial L(x)}{\partial \phi^{\alpha}(x)} -\nabla^{(x)}_{\mu}\frac{\partial L(x)}{\partial \phi^{\alpha}_{,\mu}(x)}$$ $$\tag{G}~=~\int \! d\mu(y)\left[\frac{\partial L(y)}{\partial \phi^{\alpha}(y)}\delta_M(x,y) + \frac{\partial L(y)}{\partial \phi^{\alpha}_{,\mu}(y)}d^{(y)}_{\mu}\delta_M(x,y)\right],$$ so that the 2nd functional derivative becomes symmetric $$\frac{ \delta^2_M S}{\delta\phi^{\beta}(y)\delta\phi^{\alpha}(x)} ~=~\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y) -\nabla^{(y)}_{\nu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y)\right]$$ $$+\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}_{,\mu}(y)}d^{(y)}_{\mu}\delta_M(x,y) -\nabla^{(y)}_{\nu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}_{,\mu}(y)}d^{(y)}_{\mu}\delta_M(x,y)\right]$$ $$~=~\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y) -\nabla^{(y)}_{\nu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y)\right]$$ $$-\nabla^{(x)}_{\mu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}_{,\mu}(y)}\delta_M(x,y)\right] -\nabla^{(y)}_{\nu}\nabla^{(x)}_{\mu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}_{,\mu}(y)}\delta_M(x,y)\right]$$ $$\tag{H}~=~\left[(x,\alpha) \leftrightarrow (y,\beta) \right].$$ Here we have used that $$\tag{I} [f(x)-f(y)]\delta_M(x,y)~=~0, \qquad [\nabla^{(x)}_{\mu}+d^{(y)}_{\mu}]\delta_M(x,y)~=~0.$$
References: