Finding the divergence of this vector field

Question

Let $M\in \Bbb{R}^n$ be an open set, and let $\rho: M \to \Bbb{R}$ be continuous. Then define $v: M\to \Bbb{R}^n$ as $$v(x) = \int_{r\in M} \frac{\rho(r)(x-r)}{\|x-r\|^n} \, dr$$ Where the above integral is taken component-wise. (if $n = 3$, we can imagine $\rho$ as charge density, and $v$ as the resulting electric field). I now wish to find the divergence of $v$, so first I calculate the $i^\text{th}$ partial derivative of the $i^\text{th}$ component. $$\frac{\partial v_i}{\partial x_i} = \frac{\partial}{\partial x_i} \int_{r\in M} \frac{\rho(r)(x_i-r_i)}{\|x-r\|^n} \, dr$$ Since the region of integration does not depend of $x_i$, if I'm not mistaken, we can move the derivative to the inside of the integral and differentiate using the quotient rule: $$ = \int_{r\in M}\frac{\partial}{\partial x_i}\frac{\rho(r)(x_i-r_i)}{\|x-r\|^n} \, dr = \int_{r\in M}\frac{\rho(r)\|x-r\|^n - n\rho(r)(x_i-r_i)^2\|x-r\|^{n-2}}{\|x-r\|^{2n}} \, dr$$ (since $\frac{\partial}{\partial x_i}\|x\|=\frac{x_i}{\|x\|}$) $$ = \int_{r\in M}\frac{\rho(r)\|x-r\|^2 - n\rho(r)(x_i-r_i)^2}{\|x-r\|^{n+2}} \, dr$$ So we have $$\nabla \cdot v= \int_{r\in M}\sum_{i=1}^n \frac{\rho(r)\|x-r\|^2 - n\rho(r)(x_i-r_i)^2}{\|x-r\|^{n+2}} \, dr$$ $$= \int_{r\in M} \frac{\sum_{i=1}^n(\rho(r\|x-r\|^2) - n\rho(r) \sum_{i=1}^n (x_i-r_i)^2}{\|x-r\|^{n+2}} \, dr$$ $$ = \int_{r\in M}\frac{n\rho(r)\|x-r\|^2 - n\rho(r)\|x-r\|^2}{\|x-r\|^{n+2}} \, dr = 0$$ But this seems false since it contradicts Gauss's law--we expect to get $\nabla \cdot v(x) = \lambda_{n-1}\rho(x)$ where $\lambda_{n-1}$ is the surface area of a unit $(n-1)$-sphere.

As a specific example, let $n = 1, M = (-1,1), \rho(x) = x$, and so $$v(x) = \int_{-1}^1\frac{r\cdot (x-r)}{|x-r|} \, dr = \int_{-1}^x r \, dr + \int_{x}^1 -r \, dr = x^2 - 1$$ and so $\nabla \cdot v = 2x$, which is not identically zero. Where have I gone wrong? I feel like I must be misunderstanding the conditions under which one can exchange the order of differentiation and integration, but this (first equation under the "higher dimensions" section) seems to be exactly what I'm doing (since the second term on the right hand side disappears since the area of integration does not depend on $x_i$).

Thank you in advance.

Answer 1

The problem is that the integrand is only defined if $x\ne r$. I'm not sure how this can be treated in a mathematically clean way; in physics, the usual rule is that on the differentiation the singularity gives an additional term involving the $n$-dimensional Dirac delta $\delta^{(n)}(x-r)$.

This additional delta can be nicely illustrated with the one-dimensional function (I'll use your interval, but keep $\rho$ general):

Clearly $$\frac{x-r}{|x-r|} = \begin{cases} 1 & x-r>0\\ -1 & x-r<0 \end{cases}$$ Note that, again, this function is not defined at $x=r$. However, according to the definition of the Dirac delta, $$\int_{-1}^x\delta(t-r)\,\mathrm dt = \begin{cases} 0 & x<r\ (\text{that is, } r\notin(-1,x))\\ 1 & x>r\ (\text{that is, } r\in(-1,x))\\ \end{cases}$$ and thus $$\frac{x-r}{|x-r|} = -1 + 2 \int_{-1}^x\delta(t-r)\,\mathrm dt$$ We therefore can define: $$\frac{\mathrm d}{\mathrm dx}\frac{x-r}{|x-r|} = 2\delta(x-r)$$ Thus we get \begin{aligned} \frac{\mathrm d}{\mathrm dx}v(x) &= \frac{\mathrm d}{\mathrm dx}\int_{-1}^1\frac{\rho(r)(x-r)}{|x-r|}\,\mathrm dx\\ &= \int_{-1}^1\rho(r)2\delta(x-r)\mathrm dx\\ &= 2\rho(x) \end{aligned} Inserting $\rho(x)=x$, recovers the result $\nabla v = 2x$ you calculated directly.

Answer 2

As celtschk pointed out the singularity that occurs in the integrand when $x=r$ needs to be handled carefully. Strictly speaking your manipulations were only valid for when $x\neq r$. In his book "Classical Electrodynamics" J.D Jackson presents a method for handling this singularity.

This may or may not be up to the standard of rigor you are looking for. I'm providing it in the hope that it will be useful.

---

We will "regularize" the singularity in the integrand by introducing a positive parameter $\epsilon$. In my calculations subscripts will be used to refer to the specific components of a vector; also $\partial_j=\partial/\partial x_j$.

$$v_i(x) = \int_M \left[ \frac{\rho(r)(x_i-r_i)}{(x^2+r^2-2x\cdot r + \epsilon^2)^{n/2}} \right] \, \mathrm{d}^n x $$

First we will compute $\partial_j v_i(r)$, then we will compute $\nabla\cdot v$ by contracting this against the metric $\delta_{ij}$.

$$\partial_j v_i(x) = \int_M \left[ \partial_j \frac{\rho(r)(x_i-r_i)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2}} \, \mathrm{d}^n x \right]$$

$$= \int_M\left[ \frac{\rho(r)\delta_{ij}}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2}} - \frac{(n/2)\rho(r)(x_i-r_i)(2x_j-2r_j)}{(x^2+r^2-2x\cdot r + \epsilon^2)^{n/2+1}} \right] \, \mathrm{d}^n x $$

$$= \int_M\left[ \frac{\rho(r)\delta_{ij}(x^2+r^2-2x\cdot r+\epsilon^2)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}} - \frac{(n)\rho(r)(x_i-r_i)(x_j-r_j)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}} \right]\, \mathrm{d}^n x $$

$$= \int_M\left[ \frac{\rho(r)\delta_{ij}(x^2+r^2-2x\cdot r+\epsilon^2)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}} - \frac{(n)\rho(r)(x_ix_j-r_ix_j-x_ir_j+r_ir_j)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}} \right]\, \mathrm{d}^n x $$

Now we can contract the indices, setting $i=j$ and summing from $1$ to $n$.

$$\nabla \cdot v(x) = \int_M\left[ \frac{\rho(r)n(x^2+r^2-2x\cdot r+\epsilon^2)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}} - \frac{(n)\rho(r)(x^2-2x\cdot r+r^2)}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}} \right]\, \mathrm{d}^n x $$

$$= \int_M\left[ \frac{\rho(r)\epsilon^2}{(x^2+r^2-2x\cdot r + \epsilon^2)^{n/2+1}} \right] \, \mathrm{d}^n x $$

Now when $x\neq r$ we can see that when $\epsilon \rightarrow 0$ we get a divergence of $0$ as a result. This is what we would expect in electrodynamics if $x \notin M$ because there is no charge density outside of the set. The singularity occurs when $x=r$ so lets consider the integrand in this specific case.

$$\left[ \frac{\rho(r)\epsilon^2}{(x^2+r^2-2x\cdot r + \epsilon^2)^{n/2+1}} \right]_{x=r} = \frac{\rho(x)\epsilon^2}{(\epsilon^2)^{n/2+1}}=\frac{\rho(x)}{\epsilon^n}$$

This result clearly goes to $\infty$ as $\epsilon\rightarrow 0$. Let $\epsilon_0$ be small enough that we can find an open ball in $M$ centered at $x$ for every $\epsilon < \epsilon_0$; call ball $N_\epsilon$. Then we can break the integral into the following two parts.

$$ \nabla \cdot v(x) = \int_{N_\epsilon} (\cdots) + \int_{M- N_\epsilon} (\cdots)$$

The latter part will of course go to zero as $\epsilon\rightarrow 0$. The former however will approach a finite value. For sufficiently small $\epsilon$ we will have,

$$ \int_{N_\epsilon} \left[ \frac{\rho(r)\epsilon^2}{(x^2+r^2-2x\cdot r+\epsilon^2)^{n/2+1}}\right] \mathrm{d}^n x \approx \int_{N_\epsilon} \frac{\rho(r)}{\epsilon^n} \mathrm{d}^nx = \frac{\rho(\xi)}{\epsilon^n} S_{n-1} \epsilon^n = \rho(\xi)S_{n-1} $$

where $\xi\in N_\epsilon$ and $S_{n-1}$ is the total solid angle for a $n-1$ spehere (e.g. $S_2=4\pi$). Since we have the following limit,

$$\rho(\xi)S_{n-1} \rightarrow \rho(x)S_{n-1} \text{ as } \epsilon \rightarrow 0,$$

we conclude that the sequence of integrands parameterized by $\epsilon$ defines a $\delta$ function. Informally we write,

$$ \boxed{\nabla\cdot \left[\frac{x-r}{\|x-r\|^n} \right] = S_{n-1}\delta^{(n)}(x-r)},$$

which implies that,

$$\boxed{\nabla\cdot v(x) = S_{n-1}\rho(x)}$$