There is a standard book which contains everything about electrostatics, the Laplace/Poisson equation and boundary conditions: Classical Electrodynamics by J. D. Jackson. Get the book from the library of your choice, read all chapters labeled "Electrostatics", and you will find the answers to all your questions (if you are simulating this, you need to know all this stuff anyway).
The nature of the boundary condition depends on the system you are describing. It means something else if you are calculating heat flow than electrostatics - obviously!
Let's say you have a differential equation (e.g. the Poisson equation $\Delta\varphi(\vec r)=-\frac{\varrho(\vec r)}{\varepsilon_0}$ describing electrostatics, and you solve it for the function $\varphi(\vec r)$. At the boundaries of the region (e.g. a cylinder, a cube, etc.) you have to fix some property of $\varphi(\vec r)$.
- Neumann boundary condition: You fix $\frac{\partial\varphi(\vec r)}{\partial\vec n}=\text{const}$ along the boundary, where $\vec n$ is the normal vector to the surface. It's basically the derivative of $\varphi$ if you go straight away from the surface. It can be a different value for every $\vec r$.
- Dirichlet boundary condition: You fix $\varphi(\vec r)=\text{const}$. It can be a different value for every $\vec r$.
You can only fix one of those two, or the sum (this is called Robin boundary condition).
Physical examples for electrostatics:
- Neumann boundary condition: The aforementioned derivative is constant if there is a fixed amount of charge on a surface, i.e. $\frac{\partial\varphi(\vec r)}{\partial\vec n}=\sigma(\vec r)$.
- Dirichlet boundary condition: The electrostatic potential $\varphi(\vec r)$ is fixed if you have a capacitor plate which you connected to a voltage source. E.g. if you have two capacitor plates which are at 0V and 5V, respectively, you would set $\varphi(\vec r)=0$ at the first plate and $\varphi(\vec r)=5$ at the second plate. That way, you can calculate the capacitance.
For heat flow, fixing the field $u$ (=Dirichlet BC) means fixing the temperature. If you have elements in your system that have a fixed temperature, you would use that one.
If you lenses are based on electrostatics, they probably only have Dirichlet boundary conditions, because that's how you describe a capacitor plate. If you have outer boundaries which are not capacitor plates, you should use a Neumann BC = 0 (in this case, it has nothing to do with fixing the charge), because that's the best BC for simulating an "infinite" system.
Let wave function $\Psi$ be defined on domain $D \in \mathbb{R}^n$. The Neumann condition $\frac{\partial \Psi} {\partial {\bf n}} = 0$ on the boundary $\partial D$ has a simple interpretation in terms of the probability current of $\Psi$. For $\Delta \Psi = i \partial\Psi/\partial t$ (although it's usually taken as $i \partial\Psi/\partial t = - \Delta \Psi$), the probability current at an arbitrary point ${\bf x} \in \mathbb{R}^n$ is
$$
{\bf j}({\bf x}) = i [ \Psi^*({\bf x}){\bf \nabla}\Psi({\bf x}) - \Psi({\bf x}){\bf \nabla}\Psi^*({\bf x}) ]
$$
and the normal current on $\partial D$ reads
$$
{\bf n} \cdot {\bf j} = i\; [ \Psi^* \frac{\partial \Psi}{\partial {\bf n}} - \Psi \frac{\partial \Psi^*}{\partial {\bf n}} ]
$$
(has the wrong sign, I know, but I accounted for OP's form of the Sch.eq. as $\Delta \Psi = i \partial\Psi/\partial t$).
Setting $\frac{\partial \Psi} {\partial {\bf n}} = 0$ amounts to ${\bf n} \cdot {\bf j} = 0$ everywhere on $\partial D$, thus confining the corresponding system within $D$ without an infinite potential well, as under Dirichlet conditions ($\Psi = 0$ on $\partial D$). This is the case of perfect reflection on $\partial D$.
There is a mention of this in Section 5.2 of Visual Quantum Mechanics:
Selected Topics with Computer-Generated Animations of Quantum-Mechanical Phenomena, by Bernd Thaller (Springer, 2000); Google Books link.
As for applications, one answer to another post, Can we impose a boundary condition on the derivative of the wavefunction through the physical assumptions?, pointed to the use of Neumann conditions in R-matrix scattering theory.
$$
$$
Clarification following @arivero's observation on conditions necessary to trap the system within domain D:
We can say that the system described by $\Psi$ is trapped in domain D if the total probability to locate it in D, $P_D$, is conserved in time: $dP_D/dt = 0$. In this case, if the system is initially located within D, such that $\Psi({\bf x},t=0) = 0$ for all ${\bf x} \notin D$ and $P_D(t=0) = 1$, then it will remain in D at all $t > 0$, since $P_D(t) = P_D(t=0) = 1$. If initially $P_D(t=0) < 1$ (the system has a nonzero probability to be located outside D), then we still have $P_D(t) = P_D(t=0) < 1$.
Conservation of $P_D$ is equivalent to a condition of null total probability current through the boundary $\partial D$. Note that it is not necessary to require null probability current at every point of $\partial D$, but only null total probability current through $\partial D$.
The difference can be understood in terms of path amplitudes (path-integral representation). In the former case, the amplitude $\Psi({\bf x}_1 t_1, {\bf x}_2 t_2; {\bf x} t)$ that the system "goes" from point ${\bf x}_1 \in D$ at time $t_1$ to point ${\bf x}_2 \notin D$ at time $t_2 > t_1$ while passing through point ${\bf x} \in \partial D$ at some time $t$, $t_1 < t < t_2$, is nonzero $\forall {\bf x} \in \partial D$: $\Psi({\bf x_1} t_1, {\bf x_2} t_2; {\bf x} t)≠0$. If however we demand null probability current at every point of $\partial D$, then $\Psi({\bf x}_1 t_1, {\bf x}_2 t_2; {\bf x} t)=0$, $\forall {\bf x} \in \partial D$.
In other words, null total probability current on $\partial D$ enforces weak trapping in the sense that overall $P_D(t) = $ const. and "cross-over events" across the boundary balance out. Null local probability current at every point of $\partial D$, ${\bf n} \cdot {\bf j} = 0$, corresponds to strong trapping in the sense that the system "does not cross" at any point ${\bf x} \in \partial D$. Imposing the strong trapping condition is equivalent to requiring that the weak trapping condition be satisfied by any wave function $\Psi$, as opposed to one selected $\Psi$. In this case the system is essentially confined within D at all times. Incidentally, the strong trapping condition follows from the requirement that the restriction of the system Hamiltonian on domain D remain self-adjoint.
Derivation of probability current conditions:
The free Schroedinger equation for $\Psi$, $i\partial\Psi/\partial t = \Delta\Psi$ as above (OP's choice of sign), implies local conservation of the probability density $\rho({\bf x,t}) = \Psi({\bf x},t)\Psi^*({\bf x},t)$:
$$
\frac{\partial \rho({\bf x},t)}{\partial t} + {\bf \nabla}\cdot {\bf j}({\bf x},t) = 0
$$
Integrating this over domain D yields
$$
\int_D dV\;\frac{\partial {\rho({\bf x},t)}}{\partial t} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = \frac{d}{dt}\int_DdV\;{\rho({\bf x},t)} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0
$$
which after denoting $P_D = \int_DdV\;{\rho({\bf x},t)}$ becomes
$$
\frac{dP_D}{dt} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0
$$
Imposing $dP_D/dt = 0$ necessarily means $\oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0$. Note that $\oint_{\partial D} dS\;{{\bf n}\cdot{\bf j}} = 0$ does not require ${{\bf n}\cdot{\bf j}} = 0$ at every point on $\partial D$, whereas ${{\bf n}\cdot{\bf j}} = 0$ does imply $\oint_{\partial D} dS\;{{\bf n}\cdot{\bf j}} = 0$ and $dP_D/dt = 0$.
Self-adjoint restriction of the free Hamiltonian on domain D:
A self-adjoint restriction of $H\Psi = \Delta \Psi$ on D requires that $
\int_D dV\;\Phi^* (\Delta\Psi) = \int_D dV\;(\Delta\Phi^*) \Psi$ or $
\int_D dV\;[\Phi^* (\Delta\Psi) - (\Delta\Phi)^* \Psi] = 0$.
Use $\Phi^* (\Delta\Psi) = {\bf \nabla}\cdot(\Phi^* {\bf \nabla}\Psi) - {\bf \nabla}\Phi^* \cdot {\bf \nabla}\Psi$ and Green's theorem to obtain
$$
\int_D dV\;[\Phi^* (\Delta\Psi) - (\Delta\Phi)^* \Psi] = \oint_{\partial D} dS\;[\Phi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Phi^*}{d{\bf n}}] = 0
$$
If the last condition above is to be satisfied by arbitrary $\Phi$, $\Psi$, it must hold locally:
$$
\Phi^* \frac{d \Psi}{d{\bf n}} - \Psi \frac{d \Phi^*}{d{\bf n}} = 0
$$
This means $\frac{1}{\Psi}\frac{d\Psi}{d{\bf n}} = \frac{1}{\Phi^*}\frac{d\Phi^*}{d{\bf n}} = a({\bf x}), \forall {\bf x} \in \partial D$. The case $\Phi = \Psi$ shows that $a({\bf x}) = a^*({\bf x})$. Therefore the restriction of $H$ on $D$ is self-adjoint if and only if wave functions in its support satisfy a boundary condition
$$
\frac{d\Psi}{d{\bf n}} = a({\bf x})\Psi
$$
for some given real-valued function $a({\bf x})$. In particular, this means every wave function $\Psi$ also satisfies the strong trapping condition ${\bf n} \cdot {\bf j} = 0$.
Finally note that the strong trapping condition means $\Psi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Psi^*}{d{\bf n}} = 0$, but does not imply that $\Phi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Phi^*}{d{\bf n}} = 0$ for arbitrary $\Psi$, $\Phi$. If we consider the latter expression as the matrix element of a "local current operator" ${\bf n} \cdot {\bf \hat j}$, then the strong trapping condition requires that diagonal elements of ${\bf n} \cdot {\bf \hat j}$ are 0, whereas self-adjointness requires that the entire space of wave functions is in the kernel of each ${\bf n} \cdot {\bf \hat j}$.
Best Answer
Instead of $(w, x, y, z)$, I will use $(x^a, z)$ as the $\mathbb{R}^3$ and $\mathbb{R}_+$ directions.
For the Dirichlet part, the authors are simply saying that $F_{ab} |_{z = 0} = ( \partial_a A_b - \partial_b A_a ) |_{z = 0}$ becomes zero when you plug in $A_a |_{z = 0} = \mathrm{const}$. The reverse implication would be too strong since gauge transformations can change the value of $A_a$ but not $F_{ab}$.
The Neuman part is more interesting because this requires us to go back to the action \begin{equation} S = \int_{z \geq 0} \int d^3 x -\frac{1}{4} F_{\mu\nu} F^{\mu\nu}. \end{equation} If we're careful when deriving the equations of motion, the step that involves integrating by parts produces the boundary term \begin{align} \delta S_\partial = \int d^3 x \delta A^a F_{za} \end{align} which has to vanish. If we're not using the Dirichlet condition, the only alternative is $F_{za} |_{z = 0} = 0$ which we define to be the Neuman condition.
For more about this, https://arxiv.org/abs/1902.09567 is a great modern paper.