Finding cross partial derivative of inverse function

Question

What I am interested in is to find an expression for

$$\frac{\partial^2G(S,\sigma)}{\partial S\partial\sigma}$$

where $G$ is inverse function in the first argument of function $C$ such that $S = C(G(S,\sigma), \sigma)$, and we define $C$ as

\begin{align*} C(V,\sigma) &=V N \left(d_1 \right) - D \exp \left( -rT\right) N \left(d_2\right) \\ d_1 &= \frac{\log\left(V/D\right) + \left(r + \sigma^2/2\right)T}{\sigma \sqrt{T}} \\ d_2 &= d_1 - \sigma \sqrt{T}. \end{align*}

with $S,V,T,D,\sigma\in (0,\infty)$, $r\in\mathbb{R}$, and $N$ is the standard normal cumulative distribution function. $C$ is the European call option pricing formula in the Black–Scholes model. Thus, $C: (0,\infty)^2 \rightarrow (0,\infty)$. One can show that $C$ is stricly montone in both of its arguments. There is no formula for $G$ but it can easily be computed with a bisection like method.

One can find that

\begin{align*} n(d) &= \frac{1}{\sqrt{2\sigma}}\exp(-d^2/2)\\ \frac{\partial C(V,\sigma)}{\partial \sigma} &= \sqrt{T}V n(d_1) \\ \frac{\partial^2 C(V,\sigma)}{\partial \sigma\partial V} &= \sqrt{T}n(d_1) + \sqrt{T}Vn'(d_1)\frac{\partial d_1}{\partial V} \\ &= \sqrt{T}n(d_1) + \sqrt{T}V(-d_1n(d_1)) \frac{1}{\sigma \sqrt{T}V} \\ &= \left(\sqrt{T} - d_1 / \sigma\right)n(d_1) \end{align*}

but that is not useful as far as I gather since $C$ is strictly monotone in $V$ but $\frac{\partial C(V,\sigma)}{\partial \sigma}$ is not so one cannot use results like here.

---

Ultimately, I want know if there are regions of $(S,\sigma)\in(0,\infty)^2$ where

$$ \frac{\partial^2G(S,\sigma)}{\partial S\partial\sigma} \approx 0 $$

or tend to zero. I am not sure if this is an easier question.

---

Update

One can almost get a close form solution as follows. $C$ is strictly monotone in the first argument so

$$ \frac{\partial G(S, \sigma)}{\partial S} = \left( \left.\frac{\partial C(V,\sigma)}{\partial V} \right\vert_{V = G(S,\sigma)} \right)^{-1} $$

Further, we know that

$$\frac{\partial C(V,\sigma)}{\partial V} = N(d_1(V, \sigma))$$

where we have made the depedence of $d_1$ on $(V,\sigma)$ explicit. Using the above, we find that

$$ \frac{\partial G(S, \sigma)}{\partial S} = \frac{1}{N(d_1(G(S,\sigma), \sigma))} $$

Hence,

\begin{align*} \frac{\partial G(S, \sigma)}{\partial S \partial \sigma} = - \frac{n(d_1(G(S,\sigma), \sigma))}{N(d_1(G(S,\sigma), \sigma))^2} d_1'(G(S,\sigma), \sigma) \end{align*}

where the derivative is w.r.t. $\sigma$. The only thing we are missing here is $d_1'(G(S,\sigma), \sigma)$. We can re-write $d_1$ as

$$ d_1(V,\sigma) = \frac{\log V - \log D + rT}{\sigma \sqrt{T}} + \frac{\sigma\sqrt T}{2} $$

So

$$ d_1'(G(S,\sigma), \sigma) = - \frac{\log G(S,\sigma) - \log D + rT}{\sigma^2\sqrt T} + \frac{\sqrt T}{2} + \frac{G'(S,\sigma)}{\sigma\sqrt T G(S,\sigma)} $$

This though leaves us with a $\partial G(S, \sigma) / \partial \sigma$ factor which does not appear to be nice. Finding an expression for this factor would be neat.

R code to confirm the above

#####
# assign values and functions
D     <- .8
r     <- .03
T.    <- 3

library(DtD)
G <- function(S, sigma)
  get_underlying(S = S, D = D, T. = T., r = r, vol = sigma)
Gfunc <- function(par)
  G(par["S"], par["sigma"])
d1 <- function(V, sigma)
  (log(V) - log(D) + (r + sigma^2/2) * T.) / (sigma * sqrt(T.))

S <- seq(1, 3, length.out = 20)
sigma <- seq(.1, .5, length.out = 20)
xy <- expand.grid(S = S, sigma = sigma)

#####
# compute cross partial derivative numerical differentiation
library(numDeriv)
r1 <- mapply(function(S, sigma) hessian(Gfunc, c(S = S, sigma = sigma))[2, 1],
            S = xy$S, sigma = xy$sigma)
dim(r1) <- c(length(S), length(sigma))
persp(S, sigma, r1, theta = 30, phi = 30, ticktype = "detailed")

#####
# compute cross partial derivative numerical differentiation given d G / d S
r2 <-  mapply(
  function(S, sigma)
    grad(function(sigma) pnorm(d1(G(S, sigma), sigma))^-1, sigma), 
  S = xy$S, sigma = xy$sigma)

#####
# compute cross partial derivative given almost all parts
r3 <-  mapply(
  function(S, sigma){
    V <- G(S, sigma)
    d_val <- d1(V, sigma)
    grad_fac <- grad(function(sigma) G(S, sigma), sigma)

    - (dnorm(d_val) / pnorm(d_val)^2) * (
      - (log(V) - log(D) + r * T.) / (sigma^2 * sqrt(T.)) +
        sqrt(T.) / 2
        + grad_fac / (sigma * sqrt(T.) * V))
  }, S = xy$S, sigma = xy$sigma)

#####
# check results
all.equal(c(r1), r2, tolerance = 1e-6)
#R [1] TRUE
all.equal(c(r1), r3, tolerance = 1e-6)
#R [1] TRUE