# Solving Laplace's Equation for magnetostatic problem in the absence of currents

I am attempting to simulate the resulting magnetic field distortion of paramagnetic nanoparticles in a static magnetic field in the absence of currents, in 3D. In my system, the static magnetic field is pointing in the $$\hat{z}$$ direction, so $$\vec{B}_0=\mu H_0 \hat{z}$$. I am first considering one sphere placed in the static field. I know that the solution is that of a dipole outside the sphere, and straight lines inside the sphere (as shown in Classical Electrodynamics, Jackson, Ch. 5.10). But since I am considering multiple spheres placed in different configurations, I would need to solve it numerically. I've done a bit of research and looked into a few questions on the exchange and found a few that were related but not quite what I needed:

When can Laplace equation be used to solve a magnetostatics problem

Boundary Conditions in Magnetostatics - Calculating surface current density

Magnetic field inside a perfectly diamagnetic ball

I am having trouble figuring out all the necessary boundary conditions and formulating the problem numerically. My method of attack is to numerically solve Maxwell's equations, particularly (using FreeFEM++):

\begin{align} \nabla \cdot \vec{B} &= 0\\ \nabla \times\vec{B} &= 0\\ \nabla \times\vec{H} &= 0 \end{align}

Since we are dealing with linear media, we also have:

\begin{align} \vec{H} &= \frac{1}{\mu_0}\vec{B}-\vec{M}\\ \vec{H} &= \frac{1}{\mu}\vec{B} \end{align}

Since $$\vec{H}$$ is irrotational, we can write it in terms of a magnetic scalar potential $$\phi_M$$:

$$\vec{H}=-\nabla \phi_M$$

Then, using the equations above, we can derive Poisson's equation as:

$$\nabla^2 \phi_M=\nabla \cdot \vec{M}$$

This reduces to Laplace's equation with $$\nabla \cdot \vec{M}=0$$:

$$\nabla^2 \phi_M=0$$

with the following boundary conditions:

\begin{align} \vec{H}_2 \cdot \vec{n} &= \frac{\mu_1}{\mu_2} \vec{H}_1 \cdot \vec{n}\\ \vec{H}_2 \times\vec{n} &= \vec{H}_1 \times\vec{n} \end{align}

For the potential, the first one translates to:

$$\frac{\partial}{\partial n} \phi_{M2} = \frac{\mu_1}{\mu_2} \frac{\partial}{\partial n} \phi_{M1} \tag{1}\label{eq1}\\$$

where $$\frac{\partial}{\partial n}$$ is the directional (normal) derivative.

The solver requires a variational formulation of Laplace's equation, which is given as (Jackson Ch. 1.12, using Green's Formula, with $$\delta \phi_M$$ as the test function). This functional is minimized to return $$\phi_M$$.

$$\int_V d^3r\nabla \phi_M \nabla \delta \phi_M -\int_S \delta \phi_M \frac{\partial \phi_M}{\partial n} da$$

The second term is the Neumann BC. Am I right in saying that the potential boundary condition above \eqref{eq1} is a Neumann boundary condition? There is also the issue with continuity, which is given as

$$\phi_{M1} = \phi_{M2}$$

I suppose I'm not entirely sure how to define the problem/boundaries between the sphere and the static field such that the simulation can reproduce the analytical solution. Would it be necessary to solve Laplace's equation twice (for each part) and "stitch" the fields?