# How to impose Neumann boundary conditions describing quantum Hall effect?

I am trying to solve the following differential equation : $$\nabla\cdot\left(\begin{pmatrix} c_1 & c_2 \\ -c_2 & c_1 \end{pmatrix}\nabla\phi\right) = 0$$ on a rectangular domain [1,0]x[0,1] with $$c_2$$ and $$c_1$$ constant, and $$c_2/c_1\gg 1$$.

The boundary conditions are $$\phi = \pm 1$$ on $$x=1$$ and $$x=0$$ respectively, and $$n\cdot\left(\begin{pmatrix} c_1 & c_2 \\ -c_2 & c_1 \end{pmatrix}\nabla\phi\right) = 0$$ on $$y = 0$$ and $$y=1$$.

Here's my confusion :

Since c-matrix is constant, this will evaluate to the Laplace equation, and the Neumann boundary condition will become $$\partial_y \phi = 0$$, and the information about $$c_2$$ much larger than $$c_1$$ is lost, and as expected the contour plot looks like :

sol1 = NDSolve[{D[u[x, y], x, x] + D[u[x, y], y, y] ==
NeumannValue[0, y == 0] + NeumannValue[0, y == 1],
DirichletCondition[u[x, y] == 1, x == 1],
DirichletCondition[u[x, y] == -1, x == 0]},
u, {x, y} ∈
ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}]];

ContourPlot[u[x, y] /. sol1, {x, y} ∈ mesh]


But if I set up the equation with a c-matrix then the equipotential lines look strange:

Edit: I forgot to say what omegatimestau is, it’s a function that is proportional to B, in the code snippet below it’s ~ 100, but I get the scratchy lines on the left and right sides of the box even if B = 0.

eqn = Inactive[
Div][{{1, -omegactimestau[B]}, {omegactimestau[B], 1}} .
Inactive[Grad][phi[x, y], {x, y}], {x, y}] == 0;
sol1 = NDSolve[{eqn ==
NeumannValue[0, y == 0] + NeumannValue[0, y == 1],
DirichletCondition[phi[x, y] == 1, x == 1],
DirichletCondition[phi[x, y] == -1, x == 0]},
phi, {x, y} ∈
ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}]]


The plot looks like this:

So, I have two questions :

1. How do I set up the problem so that information in the constant c-matrix is not lost?
2. What is the strange behavior happening in the second case?

• What's the definition of omegactimestau[B]? Commented Feb 26, 2022 at 2:20
• @xzczd omegatimes[B] is a function linear in B it’s a constant in the above operator definition. Sorry I didn’t specify this in the question. Commented Feb 26, 2022 at 3:02

First of all there's a simple mistake in your code, you've included == in the PDE twice! (Once when defining eqn, next inside NDSolve. ) Removing the redundant == results in a reasonable result. Still, it's slightly scratchy at the boundary, so let's make the grid denser:

<< NDSolveFEM
omegactimestau[B] = 100;
eqn = Inactive[Div][{{1, -omegactimestau[B]}, {omegactimestau[B], 1}} .
Inactive[Grad][phi[x, y], {x, y}], {x, y}] == 0;
reg = ToElementMesh[ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}],
MaxBoundaryCellMeasure -> 10^-3];

sol1 = NDSolveValue[{eqn, phi[1, y] == 1, phi[0, y] == -1},
phi, {x, y} ∈ reg];

ContourPlot[sol1[x, y], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 100]


I've omitted the zero NeumannValues, because this is the default setting of FiniteElement method.

• @xzcd Thank you for this answer, I don't understand the pdetoae function very well in the linked answer above, I have two doubts: 1) If I run the above code, it gives me error/warning ? saying 'Part 1 of {}' does not exist and 'Symbol called with 0 arguments'. 2) Do you know how I should go about solving such a differential equation because this is a well defined problem in the physical sense, by that I mean the boundary conditions have physical meaning of a constant voltage on the left and right sides of the box and current is zero on the top and bottom, and the c-matrix is conductivity. Commented Feb 26, 2022 at 21:54
• @Shubham 1. Have you executed the whole definition of pdetoae? Notice you need to execute the whole code snippet in the link. 2. If this is a well-studied problem, can you add a reference? Commented Feb 27, 2022 at 2:13
• @xzcd yes, I executed the code snippet in the link of the definition of pdetoae. The problem is well studied because the off-diagonal elements in the c-matrix are the conductivity for when voltage is applied in say, x and current is measured in y which is the Quantum Hall effect. Here's a reference, journals.aps.org/prb/abstract/10.1103/PhysRevB.47.15727. Commented Feb 27, 2022 at 2:20
• @Shubham Then I guess you don't execute definition of eqn. My eqn is the same as yours so it's omitted in the answer. Commented Feb 27, 2022 at 2:24
• @xzcd yup, it works. I think there were some conflicts in the code which gave those errors (because I had executed everything before). Anyways, thanks again for pointing out pdetoae, and the answer, as a newbie I am learning a lot. Do you have any advice for solving the pde in my question, physically it means that current goes from right to left and because the off-diagonal is large the current wants to travel along the edges. Commented Feb 27, 2022 at 2:33