How to impose Neumann boundary conditions describing quantum Hall effect?

Question

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?

Thank you for reading this.

Answer 1

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:

<< NDSolve`FEM`
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.