TITLE: Finding the Gradient from these functions.
QUESTION [1 upvotes]: Find all functions $f(x,y)$ with gradient F = $(2xy + x, x^2 + y) $
So what I did first was try to find the following equations that match the problem
$\frac{df}{dx}$ = $ (2y +1, 2x+y)$
$\frac{df}{dy}$ = $(2x+x,x^2+1)$
I did it this way because x and y are constants. From here I get lost to where to go after I take the derivatives any advice on how to proceed would be good.

REPLY [3 votes]: This problem is equivalent to finding functions $f(x,y)$ such that the differential form $\omega=(2xy+x)\mathrm dx+(x^2+y)\mathrm dy=\operatorname df$. First, we check to see if this is possible:$$
\frac \partial{\partial y}(2xy+x)=2x=\frac \partial{\partial x}(x^2+y)
$$
Thus, $\omega$ is closed and defined everywhere, so it is the differential of some function.
We can find such a function by integrating $\omega$ along the line segment from $(0,0)$ to $(x,y)$ and declaring that $f(0,0)=0$. Parametrizing this line segment as $\bar x=xt, \bar y=yt$ gives$$\begin{align}
f(x,y) &= \int_0^1[(2xyt^2+xt)x+(x^2t^2+yt)y]dt \\
&=\int_0^1[3x^2yt^2+(x^2+y^2)t]dt \\
&=x^2y+\frac12x^2+\frac12y^2.
\end{align}$$
Choosing a different point $P_a$ with $f(P_a)=0$ for the starting point of the line segment amounts to changing the constant of integration, so we have $$f(x,y)=x^2y+\frac12x^2+\frac12y^2+C.$$


TITLE: Integrating a triangle in a 2d plane, bound by $y=4x, x=1$ and $y=0$
QUESTION [1 upvotes]: I would like to complete the surface integral $\int\int x^2y^2dxdy$ bound by the triangle $y=4x$,  $y=0$,
$x=1$ integrating over y first. I can do the integral itself, however, I'm unsure of what limits to choose at each stage

REPLY [0 votes]: If you invert the integral using the Frobenius Theorem you trivially get
$$\int_{0}^{1} \int_{0}^{4x} x^2y^2dydx$$


TITLE: Help for two values of x
QUESTION [1 upvotes]: I'm looking for help. Even if you just tell me the process rather than the answer.
Given that $y=10-3x^2$, find two values of $x$ for which $y=-17$.
How would I go about answering this?

REPLY [3 votes]: Set $y = -17$, getting
$$-17 = 10 - 3x^2$$
Rearrange the equation to solve for $x^2$, and if possible take a square root.

REPLY [2 votes]: $$-17 = y = 10 - 3x^2$$
Rearranging ($x$'es to one side, the rest to the other) gives
$$3x^2 = 27$$
$$x^2 = 9$$
Taking the square root (which gives us two solutions) yields
$$x = \pm \sqrt{9} = \pm 3$$


TITLE: Need help on real analysis proof problem on sequences.
QUESTION [0 upvotes]: Prove or give a counterexample:If $ad-bc\neq0$ and $ax_n+by_n \rightarrow L$ and $cx_n+dy_n \rightarrow M$, then $x_n$ and $y_n$ both converges.
My proof:

I found a counterexample for this problem. If we let $x_n=\frac{n}{2}$
  and $y_n=-n$, then we can have $a=2, b=1$. Then, $ax_n + by_n$
  converges but $x_n$ and $y_n$ are divergent.

I am convinced my counterexample is valid, but the problem giving me $ad-bc \neq0$ makes me worried I am not doing this correctly. I feel like I have to use that inequality somehow.

REPLY [1 votes]: $$\text {We have } \quad \binom {a\;\;b}{c\;\;d} \binom {x_n}{y_n}= \binom {L+e_n}{M+f_n}$$ where $e_n \to 0$ and $f_n\to 0$ as $n\to \infty.$ Since $ad-bc\ne 0$ we have $$\binom {x_n}{y_n}=(ad-bc)^{-1} \binom {\;\;\;d\;\;\;-b}{-c\;\;\;\;\;a} \binom {L+e_n}{M+f_n}=$$ $$=(ad-bc)^{-1}\binom {\;\;\;d\;\;\;-b}{-c\;\;\;\;a}\binom {L}{M}+(ad-bc)^{-1}\binom {\;\;\;d\;\;\;-b}{-c\;\;\;\;a}\binom {e_n}{f_n}.$$


TITLE: Why is it hard to make black clothes?
QUESTION [5 upvotes]: I've noticed that we have so-called 'black' clothes everywhere, I haven't seen so far any piece of clothes that would be really black, it's always some kind of lighter black.
For instance, if you take a 'black' shirt then put some water on it, the stain will look even darker.
Why is it so hard to make (nearly) totally black clothes?

REPLY [2 votes]: This is a nontrivial problem in materials science.   People have done a lot of work both with different materials and different surface structures to create "ultra-black" absorbing structures for use in optical systems.
A big part of the problem is the statistical nature of photon absorption.  There's only a certain percent (<100) chance of absorbing a given photon with a given material.  Changing the macro structure to force a large number of reflections increases the net probability of absorption.
Oddly enough,  the military does not want perfect absorption: it's just as easy to detect a "hole" in the environment as a "beacon."  Camouflage is designed, roughly, to have a reflective parameter indistinguishable from the intended local environment.


TITLE: Distances are different by ~100-200m
QUESTION [1 upvotes]: I'm measuring distance of 2 points on Google Map and then in my program converting them into ECEF using this formula. Then using Pythagorean theorem to calculate distance between those 2 points. Difference is about 100-200m, which is fairly big for what I am doing. In my program I use double for precision.
Big question will be what method Google is using in order to calculate distance between 2 points?
Update, example:
Lon/Lat : 44.50310895983847,40.16009701404741 and 44.50151117520584,40.15672510042466. Google Earth says it's 400m.
After converting to ECEF I get 3482269.0815067333, 2938590.6893090927 and 3482536.9833229352, 2938466.2491839984. I measure distance of 276.99011m

REPLY [1 votes]: The Spherical Law of Cosines says
$$
\begin{align}
\cos\left(\frac{90^\circ}{10000\text{ km}}\text{dist}\right)
&=\sin(\text{lat}_1)\sin(\text{lat}_2)+\cos(\text{lat}_1)\cos(\text{lat}_2)\cos(\Delta\text{lon})\\
&=\cos(\Delta\text{lat})-\cos(\text{lat}_1)\cos(\text{lat}_2)(1-\cos(\Delta\text{lon}))\tag{1}
\end{align}
$$
Using the approximation $\cos(x)=1-\frac12x^2$ when $x$ is small yields
$$
\left(\frac{90^\circ}{10000\text{ km}}\text{dist}\right)^2
\stackrel.=\Delta\text{lat}^2+\cos^2(\text{lat})\Delta\text{lon}^2\tag{2}
$$
for small values of $\Delta\text{lat}$ and $\Delta\text{lon}$.
We have
$$
\begin{align}
\Delta\text{lat}&=0.0033719136227^\circ\\
\Delta\text{lon}&=0.0015977846327^\circ\\
\cos(\text{lat})\Delta\text{lon}&=0.0012211298208^\circ
\end{align}
$$
With these values, I get $\text{dist}\stackrel.=398\text{ m}$.


TITLE: Reference request: statistical inference with real world (chaotic) data, e.g. climate/weather data
QUESTION [3 upvotes]: I'm interested in the process that goes from, for example, huge data sets of temperature, pressure, precipitation etc readings to a model of the atmosphere. I'd like to know about how much relies on a priori knowledge of Navier-Stokes, and how much curve fitting is done.
Some sort of general overview of the process (ideally a fairly accessible one) would be perfect.
Edit
I'm looking for a reference, specifically. I know more or less the science from talking to people who do this stuff, but I would like to have some summary article or book chapter overview I can refer to in written work.

REPLY [3 votes]: Climate modelling is a giant science of its own, and the proportion of CFD/statistics depends on the particular model.
In general, what models do is first a simple (often uncompressable) large-scale CFD to advect the scalar fields and then apply a bunch of subrutines simulating small scale and more complex processes, like radiation transfers, heat transfers, cloud dynamics, precipitation, evatranspiration and interactions with land/ocean. Those are mainly governed by some phenomenological theories, based either on some subscale simulations, approximated theories or purely statistical models. The proper construction of that step depends on the purpose of a model and has more to do with art than science (i.e. balancing accuracy and computational time).


TITLE: Find the continuity of a function without using L hopital
QUESTION [0 upvotes]: Let [x] denote the greatest integer function & f(x) be defined in a neighbourhood of 2 by

In this we have to find the values of A & f(2) in order that f(x) may be continuous at x = 2.
I thought to use L hopital to find LHL ,but could not use it because of greatest integer function.

REPLY [1 votes]: For $1<x<2$ you have $[x+1]=2$, so the limit from the left is
$$
\lim_{x\to2^-}\frac{(\exp((x+2)\ln4))^{1/2}-16}{4^x-16}=
\lim_{x\to2^-}\frac{(\exp((x+2)\ln4))^{1/2}-16}{x-2}\frac{x-2}{4^x-16}
$$
Now,
$$
\lim_{x\to2^-}\frac{(\exp((x+2)\ln4))^{1/2}-16}{x-2}
$$
is the derivative at $2$ of the function
$$
g(x)=(\exp((x+2)\ln4))^{1/2}
$$
and
$$
g'(x)=\frac{\exp((x+2)\ln4)\ln4}{2(\exp((x+2)\ln4))^{1/2}}
$$
so
$$
g'(2)=\frac{4^4\ln4}{2(4^4)^{1/2}}=8\ln4
$$
Also
$$
\lim_{x\to2^-}\frac{4^x-16}{x-2}
$$
is the derivative at $2$ of the function $h(x)=4^x$ and $h'(x)=4^x\ln4$, so your limit is
$$
\frac{8\ln4}{16\ln4}=\frac{1}{2}
$$
For the limit from the right, do the substitution $t=x-2$, so you have
$$
\lim_{t\to0^+}A\frac{1-\cos t}{t\tan t}=
\lim_{t\to0^+}A\frac{1-\cos t}{t^2}\frac{t}{\tan t}=\dotsb
$$


TITLE: Algorithms and Simulation
QUESTION [0 upvotes]: Supposing we want to take a sample from a $N(0,1)$ distribution and i can take a sample from a $N(0,σ^2)$.
(a) Construct a disposal/rejection algorithm with function $N(0,σ^2)$, which generates a sample from the $N(0,1)$.
(b) Find the optimal value of $σ^2$.

REPLY [2 votes]: A rejection sampling algorithm to generate samples from $f$ given the ability to sample from $g$ works by generating a pair $(x,u)$ where $x\sim g$ and $u\sim U(0,1)$ and then accepting the sample iff
$$u \leq \frac{f(x)}{Mg(x)}$$
where $M$ is an upper bound on $f(x)/g(x)$. The efficiency of the algorithm is the mean number of samples accepted, which is
$$\textrm{P}\left(u \leq \frac{f(x)}{Mg(x)}\right) = \textrm{E} \left( \frac{f(x)}{Mg(x)}\right) = \int \frac{f(x)}{Mg(x)}g(x)dx = \frac{1}{M}$$
For your case, you are trying to sample from $f=N(0,1)$ by rejection sampling from $g=N(0,\sigma^2)$, which means that for $M$ you need to choose
$$M = \max_x \frac{\frac{1}{\sqrt{2\pi}}e^{-x^2/2}}{\frac{1}{\sqrt{2\pi}\sigma}e^{-x^2/(2\sigma^2)}} = \max_x \;\;\sigma e^{-\frac{1}{2}\left(1 - 1/\sigma^2\right)x^2} = \sigma$$
under the assumption that $\sigma>1$ (otherwise $f(x)/g(x)$ has no upper bound, and no rejection sampling algorithm exists). This means that the efficiency of your algorithm is
$$\frac{1}{M} = \frac{1}{\sigma}$$
which is clearly maximized when $\sigma=1$.


TITLE: Diagonalizable Matrix $A^2$
QUESTION [3 upvotes]: How can I find a matrix $A$ such that $A^2$ is diagonalizable but $A$ is not?
I have tried many different ways, but to no avail. Is there something that I am missing in the question that gives a simple answer?

REPLY [7 votes]: Take a matrix $A$ such that $A^2=0$ but $A\neq 0$.
If $A$ were diagonalizable, since the eigenvalues would be all $0$, we would have $A=0$.


TITLE: What do the two things such that "Data is fixed" and "Parameters vary" in Bayesian statistics mean?
QUESTION [6 upvotes]: While following the bayesian statistics online, the lecturer said that "Data is fixed" and "Parameters vary" in Bayesian statistics. But the explanations I got doesn't really make me understand what those things mean. The two things sound important to begin with the basic idea of Bayesian statistics. Hope to hear explanations.

REPLY [7 votes]: Under a frequentist point of view you might have an unknown parameter, say $\theta$, that you want to estimate based on some data you have collected. You assume that this true and unknown parameter is fixed. Your data are expressed through a random variable, say $X$. So, for example you are interested in maximizing a likelihood based on the probability density function $f(X\mid\theta)$. This means that you model your collected data under a belief that the probability function of your data depends on that unknown parameter. You then can estimate that parameter say by maximizing the likelihood of the data (i.e. you consider the data random, in a sense that there are a random realization of the population that you study).
Under a bayesian point of view things are a bit reversed. You do not view the parameter $\theta$ as an unknown constant, i.e. fixed at some value, that you try to estimate. You rather consider that the parameter itself has a marginal distribution $f(\theta)$ which is called a prior. This expresses you prior beliefs regarding the parameter which now is viewed as a random variable since it follows a distribution. Under such a framework you might be interested in modelling $f(\theta\mid X)$, namely update your knowledge for $\theta$ GIVEN the data you have collected. Since now the data are given, they are not random hence the "data is fixed" that your lecturer mentioned.


TITLE: Solving a system of equation modul0 5
QUESTION [2 upvotes]: Consider the system of linear equations $$\begin{pmatrix} 6 & -3\\ 2 & 6 \end{pmatrix}\begin{pmatrix} x_1\\ x_2 \end{pmatrix}=\begin{pmatrix} 3\\ 1 \end{pmatrix}  $$
a) Solve the system in $\mathbb{F}_5$
I just want to make sure my solution is correct:
We have: $$A=\begin{pmatrix} 6 & -3\\ 2 & 6 \end{pmatrix}\\ \Rightarrow A^{-1}= \frac{1}{42}\begin{pmatrix} 6 & 3\\ -2 & 6 \end{pmatrix}\\ \equiv \begin{pmatrix} 3 & \frac{3}{2}\\ -1 & 3 \end{pmatrix} \mbox{ } [5]$$
Therefore: $\begin{pmatrix} x_1\\ x_2 \end{pmatrix}=\begin{pmatrix} 3 & \frac{3}{2}\\ -1 & 3 \end{pmatrix}\begin{pmatrix} 3\\ 1 \end{pmatrix}= \begin{pmatrix} \frac{21}{2}\\ 0 \end{pmatrix} [5]   $

REPLY [4 votes]: Note that modulo $5$ one has
$$
A =
\begin{bmatrix} 6 & -3\\ 2 & 6 \end{bmatrix}
=
\begin{bmatrix} 1 & 2\\ 2 & 1 \end{bmatrix}.
$$
The latter matrix has determinant $-3 \equiv 2$, whose inverse is $3$, so
$$
A^{-1} =
\begin{bmatrix} 3 & -1\\ -1 & 3 \end{bmatrix}.
$$
Now
$$
\begin{bmatrix} x_{1}\\x_{2} \end{bmatrix}
=
A^{-1}
\begin{bmatrix} 3\\ 1 \end{bmatrix}
=
\begin{bmatrix} 3 & -1\\ -1 & 3 \end{bmatrix}.
\begin{bmatrix} 3\\ 1 \end{bmatrix}
=
\begin{bmatrix} 3\\ 0 \end{bmatrix},
$$
which is your result, as
$$
\frac{21}{2} \equiv 21 \cdot 3 \equiv 3 \pmod{5}.
$$

REPLY [1 votes]: This method is legitimate since, if a solution exists, the determinant is non-zero mod $5$, and thus is invertible mod $5$.
To continue, $21 \equiv 1 \pmod 5$ and $\frac{1}{2}=2^{-1} \equiv 3 \pmod 5$.
So the solution becomes $$\begin{bmatrix} 3 \\ 0 \\ \end{bmatrix}.$$


TITLE: Linear transformations such that $A^3=I$ but $A^2\neq I$
QUESTION [2 upvotes]: Do any such linear transformations exist? I'm new to linear algebra. Whenever I'm confronted with problems like these I tend to just guess and hope to get lucky. How would you go about constructing such a transformation?
And, to generalize, how would you construct a linear transformation such that $A^j=I$ but $A^{j-1}$ (or any smaller degree greater than one)$\neq I$?

REPLY [1 votes]: The companion matrix of the polynomial $\;x^3-1\;$ , for which this is the minimal polynomial:
$$\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\end{pmatrix}$$


TITLE: Prove $ \operatorname{im}( \phi ) \cong G/\ker( \phi)$
QUESTION [0 upvotes]: The book on abstract algebra that I'm reading uses the fact that, given the groups $G$ and $G'$ and a homomorphism $\phi:G \rightarrow G'$, then
$$ \operatorname{im}( \phi ) \cong G/\ker( \phi )$$
However the author doesn't provide a proof and simply states that it follows from "standard group theory". How could one prove the theorem above?

I would appreciate any help/thoughts!

REPLY [1 votes]: Here is a roadmap for the proof, which appears in all books:
Define $\pi : G/\ker( \phi ) \to \operatorname{im}( \phi )$ by $\pi(x \bmod \ker \phi) = \phi(x)$. Then prove:

$\pi$ is well defined
$\pi$ is a group homomorphism
$\pi$ is injective
$\pi$ is surjective


text:
\begin{document}
\maketitle

\begin{abstract}
This article presents a novel and practically useful link between geometric integration, low-discrepancy sampling and code coupling for Lagrangian and Eulerian Vlasov-Poisson solvers.\\
Low-discrepancy sequences, also called quasi-random sequences (Quasi Monte Carlo), provide convergence rates close to $\mathcal{O}( N^{-1})$ which are far superior to (pseudo) random numbers (Monte Carlo) settling in at only $\mathcal{O}(N^{-0.5})$.
Lagrangian particle methods such as PIC rely on Monte Carlo integration. The particle distributions are nonlinearly perturbed by the forward flow following the characteristics.
Hence it remains the question of whether particle methods can benefit from such quasi-random-sequences. Any nonlinear measure-preserving map does not affect the low-discrepancy of a QMC sequence such that the order of convergence remains.
This article shows that the forward flow of phase space-conserving geometric particle methods induces naturally such a measure-preserving map underlying their importance in a new framework. In this context the Hardy Krause Variation is observed to increase in the Vlasov-Poisson system for the first time. with the linear phase. Also the star discrepancy is presented for an entire PIC simulation.\\
On the other hand, Eulerian and Lagrangian solvers have different strengths and weaknesses, such that we present a novel way of transiting from a spectral discretization of the Vlasov--Poisson system to a PIC simulation. This is achieved by higher dimensional inverse transform sampling (Rosenblatt-Mück transform). In this way Markov Chain Monte Carlo techniques are circumvented which allows the use of pseudo and quasi-random numbers. In the latter case better convergence rates can be observed both in the linear and nonlinear phase. \\
\end{abstract}

Keywords: Lagrangian Particle in Cell; Vlasov--Poisson; Quasi Monte Carlo; phase space conservation; code coupling; Inverse transform sampling; Low-discrepancy; Rosenblatt-Mück transform;

\newpage

\tableofcontents

\section{Introduction}
For the solution of kinetic models such as the two dimensional Vlasov--Poisson system~\eqref{vlasov1}-\eqref{poisson1}, Eulerian solvers and Lagrangian particle methods dominate the landscape.
\begin{align}

 \partial_t f(x,v,t) + v \cdot \partial_x f(x,v,t)  + \frac{q}{m} E(x,t) \cdot \partial_v f(x,v,t)=0 \label{vlasov1}\\
 E(x,t)=- \partial_x \Phi(x,t) \\
  \partial_{xx} \Phi(x,t) = q\Big(\int_{-\infty}^{\infty} f(x,v,t)  \mathrm{d}v -1 \Big) \label{poisson1}

\end{align}
Eulerian solvers represent the phase space density $f$ on a fixed grid. Here, we consider the simplest Eulerian solver, which is based on a Fourier spectral discretization
of the entire phase space~\cite{joyce1971numerical,izrar1989integration} and a Hamiltonian splitting in time~\cite{watanabe2005vlasov}. Using the same time discretization~\cite{forest1989fourth}
but discretizing the phase space by Monte Carlo samples one can obtain the standard geometric PIC (particle in cell)\cite{LEWIS1970136,evstatiev2013variational,kraus2013variational} method.
Unfortunately, as a Monte Carlo method PIC suffers from noise~\cite{ameres2018} with the slow $N^{-\frac{1}{2}}$ convergence, which can be improved to
$N^{-(1-\epsilon) }$ by the use of Quasi Monte Carlo numbers~\cite{hickernell2005}. Still, the noise is especially a problem for small amplitudes in the initial phase of a simulation~\cite{ameres2018}.
Initially, the spectral solver appears to have no issues but after some time into the simulation, it suffers from the recurrence phenomenon \cite{einkemmer2014strategy} or
filamentations~\cite{klimas1994splitting}, which can be mitigated by anti-aliasing and Fourier filtering techniques.
Representing the density grid can be expensive in high dimensions~\cite{kormann2016sparse} and is also wasteful if large portions of phase space are practically empty. Here PIC performs better
as the markers can be placed with respect to the actual density which is known as importance sampling~\cite{caflisch98}.\\
\section{Particle in Cell}
Let us recall the fundamentals of the Particle in Cell (PIC) method~\cite{kraus2016gempic,birdsall2004plasma,decyk2011description,esirkepov2001exact}. Equation \eqref{vlasov1} describes a conservation law, which PIC solves by the methods of characteristics.
\subsection{Method of characteristics}
The characteristics $\left( \mathbf{V}(t), \mathbf{X}(t) \right)$ are defined as a curve in space time along which the value of the density $f$ stays constant:
\begin{multline}
\label{densityf:char}
 \frac{\mathrm{d}}{\mathrm{d}t} f\left(\mathbf{X}(t),\mathbf{V}(t),t\right) =
 \frac{\mathrm{d} \mathbf{X}(t) }{\mathrm{d}t}  \partial_x f \left(\mathbf{X}(t),\mathbf{V}(t),t \right)\\
 +\frac{\mathrm{d} \mathbf{V}(t) }{\mathrm{d}t}  \partial_v f \left(\mathbf{X}(t),\mathbf{V}(t),t \right)
 +\partial_t f \left(\mathbf{X}(t),\mathbf{V}(t),t \right) =0.
\end{multline}
Inserting $\partial_t f$ from~\eqref{vlasov1} into~\eqref{densityf:char} yields the equations
of motions for the characteristics of eqn.~\eqref{vlasov1}, which read
\begin{equation}
\label{vlasov1:char}
 \frac{\mathrm{d}}{\mathrm{d}t} \mathbf{V}(t) = - \frac{q}{m} E(t,\mathbf{X}(t))
 \text{ and } \frac{\mathrm{d}}{\mathrm{d}t} \mathbf{X}(t) = \mathbf{V}(t).
\end{equation}
Then $f$ as solution of eqn.~\eqref{vlasov1}
is constant along the characteristics~\eqref{vlasov1:char}, which means for given initial position in phase space
$(\mathbf{X}_0,\mathbf{V}_0)$ we have
\begin{equation}
\label{vlasov:characteristics}
 f(X(t=0),V(t=0),t=0) = f(\mathbf{X}(t),\mathbf{V}(t),t) \quad \forall t \geq 0.
\end{equation}
In this way eqn.~\eqref{vlasov1} can be solved with the method of characteristics.
Given the fields $B$ and $E$ we can follow the characteristics by solving eqn. \eqref{vlasov1:char}
with a standard ODE integrator.
We can introduce a second density $g(x,v,t)$ which solves the same Vlasov equation as $f$
\begin{equation}
\label{vlasov2}
  \partial_t g(x,v,t) + v \cdot \partial_x g(x,v,t) + \frac{q}{m}E(x,t) \cdot \partial_v g(x,v,t) =0
 \end{equation}
and call it the sampling density, prior or the law of $(\mathbf{X},\mathbf{V})$.
The initial sampling distribution $g(\cdot,\cdot,t=0)$ becomes a probability density
by imposing a normalization over the phase space $\Omega$ by $\int_{\Omega} g(x,v,t=0) \dxdv =1$ and $g(x,v,t=0)\geq 0$ for all $(x,v) \in \Omega$ .
Since $g$ follows the same Vlasov equation \eqref{vlasov2}
as $f$, see eqn.~\eqref{vlasov1}, it is constant along the same characteristics \eqref{vlasov:characteristics}.\\
The Vlasov equation $\eqref{vlasov2}$ conserves positivity and volume, therefore, $g$ stays a probability density for all $t\geq 0$, which is discussed and verified in the next section.\\
\subsection{Phase space conservation}
In order to verify that $g(x,v,t)$ is the probability density of the characteristics $(\mathbf{X}(t),\mathbf{V}(t))$ we rewrite the characteristics as a mapping $\varphi_t$.
Since $f$ is constant along the characteristics, we can implicitly define a diffeomorphism $\varphi_t: (x_0,v_0) \mapsto (x,v) $ for every $t \geq 0$ such that
\begin{align}
 \label{vlasov:characteristics:flux1}
 f(x,v,t)=f(\varphi_t(x_0,v_0),t) = f(x_0,v_0, 0).
\end{align}
The same property then also holds for $g$, namely $g(\varphi_t(x_0,v_0),t) = g(x_0,v_0, 0)$. We seek a change in variables $(x,v):=\varphi_t(x_0,v_0)$, as we are interested
in what happens with $f$ and $g$ at later times. For this denote the Jacobi determinant of $\varphi_t$ as $J_{\varphi_t}$.
In general, after a transformation has been applied onto a random deviate it's probability density has to be scaled with the according Jacobi determinant, as we will recall in the next step.
For any phase-space volume $V\subset{\Omega}$ equation~\eqref{vlasov:characteristics:flux2} then holds under the change of variables; also for $f$.
\begin{equation}
 \label{vlasov:characteristics:flux2}
 \begin{split}
 \iint_{\varphi(V)} g(x,v,t) \dxdv
 &=  \iint_{V}  g\left(\varphi_t(x_0,v_0) ,t\right)  J_{\varphi_t}(x_0,v_0) ~\mathrm{d}x_0 \mathrm{d}v_0\\
 &= \iint_{V}  g\left(x_0,v_0 ,0\right)  J_{\varphi_t}(x_0,v_0) ~\mathrm{d}x_0 \mathrm{d}v_0
 \end{split}
\end{equation}
This means that at time $t$, $(x_0,v_0)\mapsto g\left(x_0,v_0,t=0 \right) \cdot J_{\varphi_t}(x_0,v_0)$ is the probability density
for the random deviate $(\mathbf{X}(t),\mathbf{V}(t))=\varphi(X_0,V_0)$ and the Jacobian has to be taken into account.
For the Vlasov equation the Jacobi determinant is one, $J_{\varphi_t}(x,v)=1$. Hence the characteristics transport
the actual value of the probability density at every time $t$. This also holds true for a symmetric integrator, e.g.~one time step of the
symplectic Euler scheme given in equation~\eqref{vlasov:characteristics:flux:seuler}.
\begin{equation}
\label{vlasov:characteristics:flux:seuler}
\begin{split}
\varphi_t(x,v) =  \left(x + tv,  ~v + t \frac{q}{m} E(x+tv,0) \right),\\
\nabla \varphi_t(x,v)=
\begin{pmatrix}
 1 & t \\
 t \frac{q}{m} \partial_x E(x+tv,0) &  1+ t^2 \frac{q}{m} \partial_x E(x+tv,0) \\
\end{pmatrix}
\end{split}
\end{equation}
We then see that the semi-discrete flow also has the right Jacobi determinant:
\begin{equation}
\det(\nabla \varphi_t) = 1+ t^2 \frac{q}{m} \partial_x E(x+tv,0) - t^2 \frac{q}{m} \partial_x E(x+tv,0)  = 1.
\end{equation}
Yet when we consider the standard explicit Euler scheme and its Jacobi determinant
given in eqn.~\eqref{vlasov:characteristics:flux:euler} the determinant of the flow is not one.
\begin{equation}
\label{vlasov:characteristics:flux:euler}
\begin{split}
\varphi_t(x,v) &=  \left(x + tv,  ~v + t \frac{q}{m} E(x,0) \right),
\quad
\nabla \varphi_t=
\begin{pmatrix}
 1 & t \\
 t \frac{q}{m} \partial_x E(x,0) &  1\\
\end{pmatrix}
\end{split}
\end{equation}
\begin{equation}
\label{explicit:euler:jacobideterminant}
\det(\nabla \varphi_t) = 1 - t^2 \frac{q}{m} \partial_x E(x,0) \neq1
\end{equation}
Therefore the likelihood $g$ has to be rescaled accordingly such that it continuously represents the
distribution of the random deviate $(\mathbf{X}(t),\mathbf{V}(t))$. Technically, $f$ should still stay constant because the we use the method of characteristics, which leads ultimately to an inconsistency.\\
By symmetric composition it is possible to extend the dissipative explicit Euler scheme, given by $\varphi_t$ with it's left adjoint defined over the time inverted inverse $\varphi^*_t = \varphi_{-t}^{-1}$.
\begin{equation}
 \det(\nabla \varphi_t^*) = \det(\nabla \varphi_{-t}^{-1})
= \frac{1}{ \det(\nabla \varphi_{-t}) } = \frac{1}{ 1 - (-t)^2 \frac{q}{m} \partial_x E(x,0) } = \frac{1}{ \det(\nabla \varphi_{t}) }
\end{equation}
Therefore, the composition of the explicit Euler with it's adjoint should provide us a phase space conserving method. Instead of one full time step we move only a half step
and obtain a second order method. In order to work with the discrete mappings, we introduce the discrete time grid by $t_n= n\, \Delta t$. For the Vlasov--Poisson system the
electric field is obtained by the position of the particles, so it is important to note that the electric field $E(x,t_n)$ is determined from the particles $x_n$ at the $n^{\mathrm{th}}$ time step.
The discrete explicit Euler and it's adjoint read then:
\begin{align}
\label{discrete:explicit:euler}
 \varphi_{ \frac{\Delta t}{2} }&
\begin{cases}
  x_{n+1} =  x_{n+ \nicefrac{1}{2}} + \frac{\Delta t}{2} v_{n+ \nicefrac{1}{2}}\\
  v_{n+1} =  v_{n+ \nicefrac{1}{2}} + \frac{\Delta t}{2} \frac{q}{m}  E(x_{n+\nicefrac{1}{2}},t_{n+\nicefrac{1}{2}} )\\
\end{cases}\\
\label{discrete:implicit:euler}
\varphi^*_{\frac{\Delta t}{2}}= \varphi^{-1}_{-\frac{\Delta t}{2}}&

=\begin{cases}
  x_{n+ \nicefrac{1}{2}} =   x_n + \frac{\Delta t}{2} v_{n+\nicefrac{1}{2}}\\
  v_{n+ \nicefrac{1}{2}} =  v_n + \frac{\Delta t}{2} \frac{q}{m}  E( x_n + \frac{\Delta t}{2} v_{n+\nicefrac{1}{2}},t_{n+\nicefrac{1}{2}} )\\
 \end{cases}
  \end{align}
The two ways of combining \eqref{discrete:explicit:euler} and \eqref{discrete:implicit:euler} are the Crank--Nicolson,
\begin{equation}
\label{discrete:cranknicolson}
  \varphi^*_{\frac{\Delta t}{2}} \circ \varphi_{ \frac{\Delta t}{2} }
 \begin{cases}
  x_{n+1} =  x_n +  \frac{\Delta t}{2} \, \frac{v_n + v_{n+1} }{2}\\
  v_{n+1} =   v_n + \frac{\Delta t}{2} \, \frac{q}{m} \left[ E( x_n ,t_{n} ) + E( x_{n+1} ,t_{n+1} )    \right] \\
 \end{cases}
\end{equation}
and the implicit midpoint method:
\begin{equation}
\label{discrete:implicitmidpoint}
  \varphi_{ \frac{\Delta t}{2} } \circ \varphi^*_{\frac{\Delta t}{2}}
 \begin{cases}
  x_{n+ \nicefrac{1}{2}} =   x_n + \frac{\Delta t}{2} v_{n+\nicefrac{1}{2}}\\
  v_{n+ \nicefrac{1}{2}} =  v_n + \frac{\Delta t}{2} \frac{q}{m}  E( x_n + \frac{\Delta t}{2} v_{n+\nicefrac{1}{2}},t_{n+\nicefrac{1}{2}} )\\
  x_{n+1} =  x_n +  \Delta t \, v_{n+ \nicefrac{1}{2}}\\
  v_{n+1} =   v_n + \Delta t \, \frac{q}{m}  E( x_n + \frac{\Delta t}{2} v_{n+\nicefrac{1}{2}},t_{n+\nicefrac{1}{2}} )\\
 \end{cases}
\end{equation}
Both methods are fully implicit, but the crucial difference is that the implicit midpoint method approximates the electric field in the middle of the time step, whereas
the Crank--Nicolson averages the field at the beginning and the end of each time step. While both methods are considered geometric integrators with excellent long term stability
only one of them conserves the phase space volume exactly. We have already determined in eqn.~\eqref{vlasov:characteristics:flux:euler} that the explicit Euler did not conserve phase space and is dissipative,
hence it's adjoint has to be investigated. The map underlying the implicit Euler \eqref{discrete:implicit:euler} reads
\begin{equation}
 \varphi^*_t(x,v) = \left(\varphi^*_{t,x}(x,v)  ,\varphi^*_{t,v}(x,v) \right) = \left(x + t \varphi^*_{t,x}(x,v),  ~v + t \frac{q}{m} E(\varphi^*_{t,x}(x,v) ,t) \right)
 \end{equation}
Although the map is given implicitly, the Jacobi matrix can be calculated by straight forward derivation and yields also an implicit expression:
 \begin{multline}
 \nabla \varphi^*_t(x,v)=
 \begin{pmatrix}
  \partial_x \varphi^*_{t,x} & \partial_v \varphi^*_{t,x}\\
  \partial_x \varphi^*_{t,v} & \partial_v \varphi^*_{t,v}
 \end{pmatrix}\\
=
 \begin{pmatrix}
  1 + t \partial_x \varphi^*_{t,v}(x,v) & t \partial_v \varphi^*_{t,v}(x,v) \\
  t \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )  \partial_x \varphi^*_{t,x}(x,v)
  &  1 +  t \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t ) \partial_v \varphi^*_{t,x}(x,v) \\
 \end{pmatrix}
\end{multline}
By the suitable insertion of the implicitly defined derivatives the Jacobi determinant reduces to
\begin{multline}
\det(\nabla \varphi^*_t) =
    \underbrace{ \left[ 1 + t \partial_x \varphi^*_{t,v}(x,v) \right]}_{= \partial_x \varphi^*_{t,x}(x,v)}
    \left[ 1 +  t \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )  \underbrace{ \partial_v \varphi^*_{t,x}(x,v)}_{=t \partial_v \varphi^*_{t,v}(x,v)} \right]\\
  -t \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )  \partial_x \varphi^*_{t,x}(x,v) t \partial_v \varphi^*_{t,v}(x,v)= \partial_x \varphi^*_{t,x} (x,v).
\end{multline}
Unfortunately $\partial_x \varphi^*_{t,x} (x,v)$ is not known to us, but re-substituting expressions from the Jacobi matrix yields a recurrence relation
\begin{equation}
 \partial_x \varphi^*_{t,x}(x,v)=  1 + t \partial_x \varphi^*_{t,v}(x,v)  = 1+ t^2 \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )  \partial_x \varphi^*_{t,x}(x,v),
\end{equation}
which is easily resolved:
\begin{multline}
 \Rightarrow
  1 =  \partial_x \varphi^*_{t,x}(x,v) - \partial_x \varphi^*_{t,x}(x,v) t^2 \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )    \\
  \Rightarrow
  \partial_x \varphi^*_{t,x}(x,v)  = \frac{1}{1-\partial_x \varphi^*_{t,x}(x,v) t^2 \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )    }\\
  \Rightarrow
  \det(\nabla \varphi^*_t) = \frac{1}{1 - t^2 \frac{q}{m} \partial_x E( \varphi^*_{t,x}(x,v), t )}
  \label{implicit:euler:jacobideterminant}
\end{multline}
We realize that \eqref{implicit:euler:jacobideterminant} constitutes the inverse \eqref{explicit:euler:jacobideterminant} if and only if the
implicit Euler is applied first, as in \eqref{discrete:implicitmidpoint}. Note that in this case $\varphi^*_{\nicefrac{\Delta t}{2},x}(x_n,v_n) = x_{n+\nicefrac{1}{2}}$, which means
that the Jacobi determinant of the implicit midpoint scheme is one:
\begin{multline}
 \det\left[ \nabla \Big( \varphi_{\nicefrac{\Delta t}{2}} \circ \varphi^*_{ \nicefrac{\Delta t}{2}} (x_n,v_n)  \Big) \right]
 =\det\left[ \nabla \varphi_{\nicefrac{\Delta t}{2}} \big( \underbrace{\varphi^*_{ \nicefrac{\Delta t}{2}} (x_n,v_n)}_{ =( x_{n+\nicefrac{1}{2}},  v_{n+\nicefrac{1}{2}}) } \big)
  \nabla \varphi^*_{ \nicefrac{\Delta t}{2}} (x_n,v_n) \right]\\
 =   \frac{1 - \big(\frac{\Delta t}{2} \big)^2 \frac{q}{m} \partial_x E\big(  x_{n+\nicefrac{1}{2}}, t_{n+\nicefrac{1}{2}} \big)
  }{   1 - \big(\frac{\Delta t}{2} \big)^2 \frac{q}{m} \partial_x E( \varphi^*_{\nicefrac{\Delta t}{2},x}(x_n,v_n), t_{n+\nicefrac{1}{2}} )   }
  =   \frac{1 - \big(\frac{\Delta t}{2} \big)^2 \frac{q}{m} \partial_x E\big(  x_{n+\nicefrac{1}{2}}, t_{n+\nicefrac{1}{2}} \big)
  }{   1 - \big(\frac{\Delta t}{2} \big)^2 \frac{q}{m} \partial_x E(x_{n+\nicefrac{1}{2}} , t_{n+\nicefrac{1}{2}} )   } =1

\end{multline}
Contrary, for the Crank--Nicolson this holds not true. The Jacobi determinants cancel out in between half time steps, such that any series of time steps begins and ends
with a slightly dissipative half step. Hence one can call the Crank-Nicolson \textit{adjoint} phase space conserving.\\
It is important to note that the Vlasov--Poisson system is a Hamiltonian system in which our phase space coordinates $(x,v)$ coincide
with the Hamiltonian coordinates $(q,p)$. Without magnetic field the system for a single partile can be written as
\begin{equation}
 (\dot{p},\dot{q})= J^{-1} \nabla_{(p,q)} H(p,q), \quad
 J= \begin{pmatrix}
  & -I\\
  I&\\
 \end{pmatrix},
\end{equation}
with $H(p,q) = \frac{p^2}{2} + \Phi(q)$.
For different systems, we will obtain a different matrix $J$ and the coordinates $(p,q)$ cannot be identified as $(x,v)$ much longer.
An integrator is called symplectic if the mapping induced by $\varphi_t$
is symplectic with respect to $J$, which is checked by
\begin{equation}
\nabla \varphi_t(p,q)^t J \nabla \varphi_t(p,q) = J.
\end{equation}
See Hairer's lecture notes for a short introduction to Hamiltonian systems~\cite{hairer_tum2010}. Such symplectic integrators
always conserve phase-space volume and can also conserve quantities like energy but not every phase space volume-preserving integrator is symplectic, see also~\cite{mclachlan2006geometric}.
But conservation of phase space is such an important property that schemes like the Boris method perform so well although they cannot be symplectic for any system~\cite{hongqin2013}.
Many of these integrators along with detailed theory for plasma physics can already be found in~\cite{kraus2013variational}. For the Vlasov--Poisson system the commonly known schemes are symplectic Runge Kutta methods up to fourth order~\cite{forest1989fourth}, where
second order scheme corresponds to the well-known leap frog, and the first order is the symplectic Euler.
But so far it is unclear whether symplecticity provides advantages concerning the conservation of low-discrepancy, which is why we restrict ourselves to phase space conservation.

\subsection{Monte Carlo integration with particles}
So far we did not address how to solve the actual Poisson equation and acquire the electric fields.
We will now slightly deviate in notation from the standard Particle-In-Cell (PIC) method~\cite{birdsall2004plasma}. The introduction of the probability density function $g$ allows us to define the characteristics $\mathbf{X}(t)$ and $\mathbf{V}(t)$ as random variables for each time $t$, such that the
trajectories in time form a stochastic process~\cite{oksendal2003stochastic}
describing the solution to eqn.~\eqref{vlasov1}.\\
Before solving the Poisson equation one can explain the stochastic setting by estimating the kinetic energy $\mathcal{H}_T$ which is a moment of the solution $f$. The characteristics are random variables with joint probability density $g$. Suppose $g$ supports $f$ which means $ \mathrm{supp}(f) \subset \mathrm{supp}(g) $. Then by inserting $g$ integrals over $f$ can be linked to expected values over $X$ and $V$.
\begin{equation}
\label{estimator:ekin1}
\begin{split}
\mathcal{H}_T(t) &= \frac{1}{2} \int_{\Omega} v^2 \, f(x,v,t) \dxdv \\
& = \frac{1}{2} \int_{\Omega} v^2 \, \frac{f(x,v,t)}{g(x,v,t)}~g(x,v,t) \dxdv \\
 &=\frac{1}{2} \,\mathbb{E}\left[ \mathbf{V}(t)^2 \frac{f(\mathbf{X}(t),\mathbf{V}(t),t)}{g(\mathbf{X}(t),\mathbf{V}(t),t)} \right] \\
 \end{split}
\end{equation}
The values of $g$ and $f$ over time are constant along the characteristics (see eqn.~\eqref{vlasov:characteristics})
such that eqn.~\eqref{estimator:ekin1} simplifies to eqn.~\eqref{estimator:ekin2}.
\begin{equation}
\label{estimator:ekin2}
\mathcal{H}_T(t) =
\frac{1}{2} \,\mathbb{E}\left[ \mathbf{V}(t)^2 \frac{f(\mathbf{X}(t),\mathbf{V}(t),t)}{g(\mathbf{X}(t),\mathbf{V}(t),t)} \right]
=\frac{1}{2}\, \mathbb{E}\left[ \mathbf{V}(t)^2 \frac{f(X(0),V(0),0)}{g(X(0),V(0),0)} \right] \\
\end{equation}
In order to get an estimate of the expectation in eqn.~\eqref{estimator:ekin2} one has to use Monte Carlo integration. We define $N_p$ independently and identically distributed (i.i.d.) samples
$\left(\mathbf{x}^0_k, \mathbf{v}^0_k \right)_{k=1,\dots,N_p}$ of the random deviates $(X(0),V(0))$ using the knowledge of the probability density $g(x,v,t=0)$. These samples are called markers or particles.
The samples $\left(\mathbf{x}^t_k, \mathbf{v}^t_k \right)_{k=1,\dots,N_p}$ can be advanced over time using a suitable phase space conserving time integrator. Then at any point in time they are distributed according to $g(\cdot,\cdot, t)$ as a solution to the Vlasov equation.
Note that the plasma likelihood
\begin{equation}
\mathbf{f}^t_k := f( \mathbf{x}_k^t, \mathbf{v}_k^t,t) = f( \mathbf{x}_k^0, \mathbf{v}_k^0,0)      =\mathbf{f}^0_k
\end{equation}
and the sampling likelihood
\begin{equation}
\mathbf{g}^t_k := g( \mathbf{x}_k^t, \mathbf{v}_k^t,t) = g( \mathbf{x}_k^0, \mathbf{v}_k^0,0) \mathbf{g}^0_k
\end{equation}
stay constant over time. In the common notation of collisionless PIC schemes
the ratio between those two likelihoods is referred to as the time-independent particle weight $\mathbf{w}_k^0$.
\begin{equation}
 \mathbf{w}_k= \mathbf{w}_k^t = \frac{\mathbf{f}^t_k }{\mathbf{g}^t_k } = \frac{\mathbf{f}^0_k }{\mathbf{g}^0_k } =  \mathbf{w}_k^0
\end{equation}
This allows us to estimate the kinetic energy using the samples in eqn.~\eqref{estimator:ekin3}.
\begin{equation}
\label{estimator:ekin3}
\begin{split}
\mathcal{H}_T(t) &=\frac{1}{2}\, \mathbb{E}\left[ \mathbf{V}(t)^2 \frac{f(X(0),V(0),0)}{g(X(0),V(0),0)} \right]\\
&\approx \frac{1}{2}\,\frac{1}{N_p} \sum_{k=1}^{N_p} \big(\mathbf{v}_k^t\big)^2 \frac{\mathbf{f}^0_k }{\mathbf{g}^0_k }
= \frac{1}{2}\, \frac{1}{N_p} \sum_{k=1}^{N_p} \big(\mathbf{v}_k^t\big)^2 \mathbf{w}_k.
\end{split}
\end{equation}
Note that sometimes the factor $\frac{1}{N_p}$ is included into the particle weight $\mathbf{w}_k$.
\subsection{Weak Poisson solve with particles}
It only remains to solve the Poisson equation~\eqref{poisson1} using the samples $\left(\mathbf{x}^0_k, \mathbf{v}^0_k \right)_{k=1,\dots,N_p}$. This is commonly done in weak form given a test function $\text{ for } \varphi \in H_1([0,L])$, see also \cite{kraus2016gempic}.
\begin{equation}
 - \int_0^L \Phi(x,t)  \varphi(x) \dx
 = \int_0^L q\Big(\int_{-\infty}^{\infty} f(x,v,t)   \mathrm{d}v -1 \Big) \varphi(x) \dx
\end{equation}
The only unknown is the right hand side,
\begin{equation}
q \int_0^L \int_{-\infty}^{\infty} f(x,v,t)   \varphi(x) \,\mathrm{d}v \mathrm{d} x
\end{equation}
which depends on the density $f$. But in the same manner as we calculated the
kinetic energy we can use Monte Carlo estimator to calculate the right hand side.
\begin{equation}
\begin{split}
q \int_0^L \int_{-\infty}^{\infty} f(x,v,t)   \varphi(x) \,\mathrm{d}v \mathrm{d} x
&=  q\,\mathbb{E} \left[  \varphi\big(\mathbf{X}(t)\big) \frac{f(\mathbf{X}(t),\mathbf{V}(t),t)}{g(\mathbf{X}(t),\mathbf{V}(t),t)} \right]\\
&\approx
 q\,\frac{1}{N_p} \sum_{k=1}^{N_p} \big(\mathbf{v}_k^t\big)^2
\underbrace{\frac{\mathbf{f}^0_k }{\mathbf{g}^0_k } }_{=\mathbf{w}^0_k}
\end{split}
\end{equation}
In order to obtain the common Particle-in-Cell formulation a spline basis as in eqn.~\eqref{linearsplinebasis} can be used for the test function $\varphi$ and the solution $\Phi$.
The electric field is then easily obtained $E(x,t)= - \nabla \Phi(x,t)$ at any time such that it can be used in the phase-space conserving ODE integrators for advancing the particles.

\section{Monte Carlo and Quasi Monte Carlo}
As commonly known the Monte Carlo estimator converges with $\mathcal{O}(N^{-0.5})$ when using $N$ pseudo random samples. Samples obtained from low discrepancy sequences,
also known as Quasi Monte Carlo (QMC), can obtain convergence up to $\mathcal{O}(N^{-1} \log(N)^{d-1}  )$ under certain smoothness conditions onto the integrand~\cite{morokoff1994quasi,basu2016transformations}. Also, the convergence rate is not restricted to the Lesbegue measure~\cite{aistleitner2016tusnady}.
Yet the measure of error for integration with these low discrepancy sequences is the Hardy-Krause variation, see~\cite{hickernell2005} for an overview and also \cite{drmota2006sequences,hlawka1971definition}.
High order scrambling by Dick, see~\cite{dick2011higher,dick2013high}, leads to convergence rates up to $\frac{7}{2}$ but requires even smoother integrands.

\begin{equation}
\label{eq:Koksma-Hlawka}
 \Big| \int f(x) \mathrm{d}x - \frac{1}{n} \sum_{p=1}^n f(x_p)   \Big| \leq V(f) D_n^*
\end{equation}
where $V(f)$ denotes the Hardy Krause variation of $f$ and $D_n^*$ the star discrepancy of the point set $x_n$.
The total variation in the sense of Hardy Krause depends on the scale of the smallest features of $f$ and can be calculated for continuously differentiable functions, see \cite{morokoff1994quasi}.
Almost any randomly chosen sequence will be bounded as
\begin{equation}
  D_n = \mathcal{O} \Bigg( \sqrt{ \frac{ \log( \log(n ))}{n}} \Bigg) \text{ with } D_n^* \leq  D_n \leq 2^d D_n^*,
\end{equation}
hence the observed $\mathcal{O}\sqrt{n}$ convergence. On the other hand the Halton and Sobol quasi random sequences yield a asymptotically smaller star discrepancy:
\begin{equation}
 D_n^* = C_d \frac{\log(n)^d}{n} + \mathcal{O}\Bigg( \frac{ \log(n )^{d-1}}{n}  \Bigg) ,  \quad C_d \ge 0.
\end{equation}
For an overview over different discrepancies and sequences~\cite{morokoff1994quasi} is recommended. For the Sobol sequence the best upper bound
for $C_d$ is $5.28$. But in~\cite{thiemard2001optimal} it is shown that the current theory on those bounds is insufficient for any practical use. Note that the total variation is defined as
\begin{equation}
 V(f) = \int | \nabla f(x)  | \mathrm{d} x,
\end{equation}
but in the world of low discrepancy sequences one mostly uses the definition in the sense of Hardy and Krause, which is the sum of the $L^1$ norm of all first order partial derivatives, except that
the identically mixed derivatives are only accounted once. For our two dimensional phase space this reads
\begin{equation}

\label{tvhk2d}
  V(f) =  \iint | \partial_{x} f(x,v,t) | \, \mathrm{d}v\mathrm{d}x
  \,+\,  \iint | \partial_{v} f(x,v,t) | \, \mathrm{d}v\mathrm{d}x
  \,+\,  \iint | \partial_{x} \partial_{v} f(x,v,t) | \, \mathrm{d}v\mathrm{d}x.
\end{equation}
Naturally, the question arises how the total variation changes over time in the Vlasov--Poisson system. Since it is quite hard to compute~\eqref{tvhk2d} from a marker distribution we used a
pseudo spectral solver. Here fig.~\ref{fig:tv_hk_overview} clearly shows increases in the variation of several orders of magnitudes due to the development of small scales. This means
according to the Koksma-Hlawka inequality~\ref{eq:Koksma-Hlawka}, that even if the discrepancy of the markers stays constant much more markers are needed in the nonlinear phase.
\begin{figure}[H]
\centering
 \includegraphics[width=0.8\textwidth]{./plots/hardykrause_variation/tv_hk_overview.eps}
 \caption{Total variation of the distribution function $f$ in the sense of Hardy and Krause for different test cases of the Vlasov Poisson. Using
           pseudo spectral solver with a third order symplectic Runge Kutta scheme at $\Delta t = 0.05$ and high resolution $N_x=512, ~N_v=1024$ the total variation
           is calculated at each time step according to eqn.~\eqref{tvhk2d}. For linear Landau damping, there are only minor changes in the distribution function since only the recurrence phenomenon is causing a
           long time disturbance. Both nonlinear Landau damping and the Bump-on-tail instability create many small scale features in the distribution function, which leads to the harsh increase in variation.
           The vortex in the Bump-on-tail instability, also known as the BGK mode, take some time to fully develop whereas the nonlinear Landau damping ``folds''
           the phase space much quicker, which explains the difference in the initial development of the variation.}
 \label{fig:tv_hk_overview}
\end{figure}
Therefore also the star discrepancy of the markers in a PIC simulation has to be investigated more closely.
There are various ways of estimating the (star) discrepancy of a given point set~\cite{doerr2014calculation}, where we chose the method provided by~\cite{thiemard2001optimal}.
The original implementation provided by~\cite{thiemard2001optimal} is restricted to computing the discrepancy in a quadratic box. Hence we periodically estimate the
star discrepancy of all markers in the phase space box $ (x,v) \in [0,2]\times[-1,1] $. Originally the markers are sampled uniformly in phase space, such that we expect this
uniformity to be preserved over time in any subinterval. This way, unfortunately, the sampling has to be cut off in the velocity domain such that we chose $|v|\leq 8$ initially.
The testing box is chosen smaller to avoid the influence of this boundary since it is not guaranteed that the support of the sampling distribution $g$ stays constant.
We already know that the standard Euler is dissipative, so we expect it to influence the star discrepancy. As can be seen fig.~\ref{fig:pic:star_discrepancy}
the symplectic Euler suffers only from a minor change in the discrepancy whereas the standard Euler suffers from an increase of orders of magnitude.
\begin{figure}[H]
\centering
\begin{tabular}{c | c}
 \hline
 \multicolumn{2}{c}{Pseudo random}\\
 \includegraphics[width=0.48\textwidth]{./plots/star_discrepancy/landau_strong_rnd_star_discrepancy.png}
 &
 \includegraphics[width=0.48\textwidth]{./plots/star_discrepancy/bumpontail_rnd_star_discrepancy.png}\\
 \hline
  \multicolumn{2}{c}{Quasi Monte Carlo (Sobol)}\\
 \includegraphics[width=0.48\textwidth]{./plots/star_discrepancy/landau_strong_qmc_star_discrepancy.png}&
 \includegraphics[width=0.48\textwidth]{./plots/star_discrepancy/bumpontail_qmc_star_discrepancy.png}\\
  \hline
 (a) nonlinear Landau damping  & (b) Bump-on-tail instability\\
 \end{tabular}

\caption{The star discrepancy in a long term PIC simulation of nonlinear Landau damping (a) and a Bump-on-tail instability (b) with a total of $N_p=2\cdot10^6$ particles sampled uniformly from the Sobol sequence. The explicit euler ``euler''
          is dissipative which worsens the discrepancy compared to the phase space conserving symplectic euler ``seuler''.}
\label{fig:pic:star_discrepancy}
\end{figure}
In view of Koksma-Hlawka inequality~\eqref{eq:Koksma-Hlawka} these results strongly recommend the use of uniformity preserving methods, such as the symplectic Euler because otherwise much more markers are required.

\subsection{Inverse Transform Sampling}
For a given phase space density $f(x,v)$ we can define the sampling density $g$ as
\begin{equation}
\label{normalizef}
 g_{\mathbf{X},\mathbf{V}}(x,v)= \frac{|f(x,v)|}{ \int_{ v_{\min}}^{v_{\max}} \int_{x_{\min}}^{x_{\max}} |f(x,v)|~  \mathrm{d}x\mathrm{d}v }
\end{equation}
There are various ways of sampling from an arbitrary probability density $g$. For pseudo random numbers popular choices are Markov Chain Monte Carlo, Gibbs sampling or even
the inefficient rejection sampling~\cite{caflisch98,liu2008monte,owen2013,christian2007monte}.
Although there exist MCMC algorithms~\cite{owen2005quasi}, the Monte Carlo schemes do not easily extend to
low discrepancy sequences.\\
Caflisch~\cite{caflisch98} already notes that the simplest way of sampling from both pseudo- and quasi-random numbers is inverse transform sampling (ITS) using the
inverse cumulative probability density. This method, especially for higher dimensions is also known as Rosenblatt-Mück transformation~\cite{rosenblatt1952remarks,Hlawka1972}.\\
For a one dimensional probability density $p_{\mathbf{X}} : [x_{\min},x_{\max}] \rightarrow \mathbb{R}$ the corresponding cumulative density reads
\begin{equation}
 P_{\mathbf{X}} (x) = \int_{x_{\min}}^{x}  p(x) \mathrm{d} x.
\end{equation}
Given a uniformly distributed pseudo- or quasi-random number $\mathbf{u} \sim \mathcal{U}(0,1)$ the corresponding sample $\mathbf{x}$
from the probability density $p_{\mathbf{X}}$ is obtained by using the inverse cumulative distribution function $ P^{-1}_{\mathbf{X}} $
\begin{equation}
 P_{\mathbf{X}} (\mathbf{x}) = \mathbf{u}  \quad \Rightarrow \quad \mathbf{x} =  P^{-1}_{\mathbf{X}} (\mathbf{u}).
\end{equation}
Note that it is also common to solve the inversion by Picard iterations or a Newton method.\\
For a given initial condition $f(x,v,t=0)$ to be used in a PIC code the density is mostly so simple that
it decomposes into a tensor product of one dimensional pieces~\cite{caflisch98}, which are then sampled by the one dimensional ITS. \\
We know that PIC codes perform poorly in situations with small perturbations, which is mostly the case in the initial phase of a simulations. Therefore, one might
use a spectral solver to start the simulation for $t \in [0, t_0]$  and then continue with PIC for $t \in [t_0, t_{\max}]$ by using the density
$f(x,v,t=t_0)$ as the initial condition. This requires importance sampling from the density $f(x,v,t=t_0)$
which cannot be done anymore by ITS one dimensional pieces. Therefore, we introduce two dimensional (inverse) transform sampling~\cite{olver2013fast} for sampling from
an arbitrary PDF $ g_{\mathbf{X},\mathbf{V}}(x,v)$.\\
On starts with a sample in the first dimension, which is obtained by ITS from a marginal density. Then as we walk through the dimensions we inverse transform sample
from the conditional marginal distribution (integrating over all the higher dimension)
\textit{given} all the previous samples from the lower dimensions.\\
In the two dimensional case we are given a uniformly distributed pseudo- or quasi-random tuple $(\mathbf{u}_x,\mathbf{u}_v) \in \mathcal{U} (0,1)^2$ and have to obtain
the sample $(\mathbf{x},\mathbf{v})$. The first marginal distribution reads
\begin{equation}
 g_{\mathbf{X}} (x) = \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(x,v)~  \mathrm{d}v,
\end{equation}
and can be sampled from by finding $ \mathbf{x}$ such that $G_{\mathbf{X}} (\mathbf{x})=\mathbf{u}_{x}$, where
\begin{equation}
  G_{\mathbf{X}} (x) =  \int_{ v_{\min}}^{v_{\max}} \int_0^x  g_{\mathbf{X},\mathbf{V}}( \hat{x},v)~  \mathrm{d}\hat{x}\mathrm{d}v
\end{equation}
This is nothing else than using $\mathbf{u}_x$ for ITS from the marginal density $g_{\mathbf{X}}$.
Now \textit{given} the sample $\mathbf{x}$ the conditional density for the likelihood of having a particle at $v$ reads
\begin{equation}
 g_{\mathbf{X}= \mathbf{x}, \mathbf{V}} (v) = \frac{  g_{\mathbf{X},\mathbf{V}}(\mathbf{x},v)  }{ g_{\mathbf{X}} (\mathbf{x})  }
 =
 \frac{  g_{\mathbf{X},\mathbf{V}}(\mathbf{x},v)  }{ \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(\mathbf{x},v)~  \mathrm{d}v  }
\end{equation}
Note that $g_{\mathbf{X},\mathbf{V}}$ gives as the probability of having both $(x,v)$, but since we already fixed $x=\mathbf{x}$ we have to normalize
with the corresponding probability for $\mathbf{x}$, namely $g_{\mathbf{X}}(\mathbf{x})$.
After inverting the corresponding cumulative conditional probability density,
\begin{equation}
  G_{\mathbf{X}= \mathbf{x}, \mathbf{V}} (v) =  \int_{ v_{\min}}^{v}  ~ g_{\mathbf{X}= \mathbf{x}, \mathbf{V}} ( \hat{v})~\mathrm{d}\hat{v}
  =
    \frac{  \int_{ v_{\min}}^{v}  g_{\mathbf{X},\mathbf{V}}(\mathbf{x},\hat{v})  ~\mathrm{d}\hat{v} }{ \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(\mathbf{x},\hat{v})~  \mathrm{d}\hat{v}  }

\end{equation}
according to
\begin{equation}
G_{\mathbf{X}= \mathbf{x}, \mathbf{V}} (\mathbf{v}) = \mathbf{u}_v
\end{equation}
the second sample $\mathbf{v}$ is obtained. Considering the following map
\begin{equation}
 \begin{split}
\Pi:~& [x_{\min},x_{\max}] \times [v_{\min},v_{\max}] \rightarrow [0,1]^2\\
  &(x,v) \mapsto
 \begin{pmatrix}
   G_\mathbf{X}( x) \\     G_{\mathbf{X}= x, \mathbf{V}} (v )
 \end{pmatrix}
 =
\begin{pmatrix}
    \int_{ v_{\min}}^{v_{\max}} \int_0^x  g_{\mathbf{X},\mathbf{V}}( \hat{x},\hat{v})~  \mathrm{d}\hat{x}\mathrm{d} \hat{v} \\
    \frac{  \int_{ v_{\min}}^{v}  g_{\mathbf{X},\mathbf{V}}(x,\hat{v})  ~\mathrm{d}\hat{v} }{ \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(x,\hat{v})~  \mathrm{d}\hat{v}  }.
 \end{pmatrix}
 \end{split}
\end{equation}
Then the entire procedure of inverse transform sampling can be described by the inverse map
$\Pi^{-1}:~  [0,1]^2 \rightarrow [x_{\min},x_{\max}] \times [v_{\min},v_{\max}]$. By considering the Jacobi matrix of $\Pi$,
\begin{multline}
\mathrm{D} \Pi (x,v)=\\
 \begin{pmatrix}
    \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}( x,\hat{v})\mathrm{d} \hat{v}      & 0 \\
    \frac{
\int_{ v_{\min}}^{v}  \partial_x  g_{\mathbf{X},\mathbf{V}}(x ,\hat{v})  \mathrm{d}\hat{v}
\int_{ v_{\min}}^{v_{\max}} g_{\mathbf{X},\mathbf{V}}(x,\hat{v})  \mathrm{d}\hat{v}
   -
\int_{ v_{\min}}^{v}  g_{\mathbf{X},\mathbf{V}}(x,\hat{v})  \mathrm{d}\hat{v}
\int_{ v_{\min}}^{v_{\max}}   \partial_xg_{\mathbf{X},\mathbf{V}}(x,\hat{v})  \mathrm{d}\hat{v}
}{
    \left( \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(x,\hat{v})  \mathrm{d}\hat{v} \right)^2
    }
    &\frac{  g_{\mathbf{X},\mathbf{V}}(x,v) }{ \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(x,\hat{v})  \mathrm{d}\hat{v}  }
 \end{pmatrix}
\end{multline}
the Jacobian introduced by the map $\Pi^{-1}$ is then the Jacobi determinant of $\Pi$:
\begin{equation}
 \mathrm{det} \left(  \mathrm{D} \Pi(x,v)  \right)  =
 \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}( x,\hat{v})~\mathrm{d} \hat{v}
 \frac{  g_{\mathbf{X},\mathbf{V}}(x,v) }{ \int_{ v_{\min}}^{v_{\max}}  g_{\mathbf{X},\mathbf{V}}(x,\hat{v})~  \mathrm{d}\hat{v}}
 =g_{\mathbf{X},\mathbf{V}}(x,v).
\end{equation}
This proofs that transforming uniform samples by $\Pi^{-1}$ introduces the Jacobian $g_{\mathbf{X},\mathbf{V}}$ which means we sample from $g_{\mathbf{X},\mathbf{V}}$, which is what we wanted.
The discrete map, based on bilinear interpolation, satisfies smoothness conditions such that the sampling is valid also for QMC numbers, see also~\cite{aistleitner2013low,basu2016transformations}.\\
Suppose a probability density is given in spectral form as
\begin{equation}
\label{fourierg}
  g(x,v) = \sum_{k_x}  \sum_{k_v} \hat{g}(k_x,k_v) \euler^{\imagi ( k_x  (x-x_{\min}) + k_v (v-v_{\min}) )},
\end{equation}
then the marginals can be represented directly by
 \begin{equation}
  G_{\mathbf{X}} (x)  = \hat{g}(0,0) x +   \sum_{k_x\neq 0}  \hat{g}(k_x,0) \frac{1}{ \imagi k_x} \left( \euler^{\imagi  k_x  x-x_{\min})  } -1 \right).
 \end{equation}
Apart from the fact that such dense Fourier interpolation is very expensive another complication arises in the combination with the used spectral solvers.
The obtained distribution function $f$ is not necessarily non-negative such that by the normalization in eqn.~\eqref{normalizef} an absolute value is introduced
in eqn.~\eqref{fourierg} which makes the corresponding anti-derivatives much more complicated. Although it might be physically not reasonable particles can still have the negative weight
$\mathbf{w} = \frac{f( \mathbf{x}, \mathbf{v} ) }{g_{\mathbf{X},\mathbf{V}(\mathbf{x}, \mathbf{v} ) }}$ in order to continue the PIC as close as possible to the spectral solution.
Also, anti-aliasing and Fourier filtering can mitigate the problem but there are no guarantees.
But be warned, ignoring the non-negativity constraint for the sampling density leads to non-monotonic increasing marginal densities
such that the inversion is not well-posed anymore such that this will not work.
Since the Fourier interpolation is expensive anyway and the density $f$ is with the help of the FFT available on a Cartesian grid anyhow, interpolation suggests itself.
\subsection{Inverse Transform Sampling from a Bilinear Interpolant}

Suppose a two dimensional sampling density $g$ is given as a bilinear interpolant at grid points $(x_i, y_j)$ with values $\left( g_{\mathbf{X,Y}}^{i,j} \right) = g_{\mathbf{X,Y}}(x_i,y_j) $ as
\begin{equation}
\begin{split}
g_{\mathbf{X,Y}}(x,y)=
\begin{bmatrix}
 1- \frac{ x - x_{i,j}}{\Delta x}  & \frac{ x - x_{i,j}}{\Delta x}
\end{bmatrix}
\begin{bmatrix}
 g^{i,j} & g^{i,j+1}\\
 g^{i+1,j} & g^{i+1,j+1}\\
\end{bmatrix}
\begin{bmatrix}
 1- \frac{ y - y_{i,j}}{\Delta y}  & \frac{ y - y_{i,j}}{\Delta y}
\end{bmatrix}\\
\text{ for } x \in [ x_i, x_{i+1}] \text{ and } y \in [ y_j, y_{j+1}]

\end{split}.
\end{equation}
Then the integral over the entire domain is given by the trapezoidal rule
\begin{equation}
  \int_{y_{\min}}^{y_{\max}} \int_{x_{\min}}^{x_{\max}}  g_{\mathbf{X,Y}}(x,y) ~\mathrm{d}x \mathrm{d}y =
  \sum_{i=1}^{N_x-1}  \sum_{j=1}^{N_y-1} \left(g^{i,j} + g^{i+1,j} +g^{i,j+1} +g^{i+1,j+1} \right) \frac{\Delta x \Delta y}{4} = 1.
\end{equation}
In the case that $(g_{\mathbf{X,Y}}^{i,j})\geq 0$ for all $i,j$ the bilinear interpolation guarantees the positivity of the interpolant $g_{\mathbf{X,Y}}$.
Provided a pair of uniform random or quasi-random numbers $ (u_x,u_y) \in [0,1]^2$ we describe the step by step procedure for determining the unique sample $(\mathbf{x},\mathbf{y})$
that corresponds to inverse transform sampling of the density according by inverse transform sampling of the density $g_{\mathbf{X,Y}}$.
For sampling in the first dimension we need the marginal density
\begin{equation}
 g_{\mathbf{X}} (x)  = \begin{bmatrix}
 1- \frac{ x - x_{i}}{\Delta x}  & \frac{ x - x_{i}}{\Delta x}
\end{bmatrix}
\begin{bmatrix}
g_{\mathbf{X}}^i \\ g_{\mathbf{X}}^{i+1}
\end{bmatrix}
 \text{ for } x \in [ x_i, x_{i+1}]
\end{equation}
where the $g_{\mathbf{X}}$ at the grid points is exactly given by the trapezoidal rule as
\begin{equation}
g_{\mathbf{X}}^i:= g_{\mathbf{X}}(x_i) =  \int_{y_{\min}}^{y_{\max}}   g_{\mathbf{X,Y}}(x_i,y) ~\mathrm{d}y
=   \sum_{j=1}^{N_y-1} \frac{  \left( g^{i,j+1} +g^{i,j+1} \right)\Delta y}{2}.
\end{equation}
This allows us to calculate the cumulative distribution function of the marginal density in $\mathbf{X}$ as
\begin{equation}
  G_{\mathbf{X}} (x) =  \Delta x
  \begin{bmatrix}
\frac{ x - x_{i}}{\Delta x} (  1- \frac{ x - x_{i}}{2 \Delta x} )  &  \left(\frac{ x - x_{i}}{\Delta x} \right)^2
\end{bmatrix}
\begin{bmatrix}
g_{\mathbf{X}}^i \\ g_{\mathbf{X}}^{i+1}
\end{bmatrix}
+ \sum_{1 \leq k<i } g_{\mathbf{X}}^i
\end{equation}
Given a uniform random number $ \mathbf{u}_x \in [0,1]$, we search for $\mathbf{x}$ such that
\begin{equation}
 G_{\mathbf{X}} ( \mathbf{x} ) = \mathbf{u}_x.
\end{equation}
Since the underlying polynomial is only of quadratic type and monotonic increasing this inversion can be solved directly.
The cell index $i$ of $\mathbf{x}$ is found by determining the largest $i$ such that
\begin{equation}
 \sum_{1 \leq k<i } g_{\mathbf{X}}^i  < \mathbf{u}_x,
\end{equation}
which gives
\begin{equation}
\label{eqn:bilinear:cellinversionx}
 \mathbf{x} = x_i + \Delta x
 \frac{ - g_{\mathbf{X}}^i + \sqrt{ ( g_{\mathbf{X}}^i )^2 + 2 (g_{\mathbf{X}}^{i+1} - g_{\mathbf{X}}^i )
\left( \mathbf{u}_x - \sum_{1 \leq k<i } g_{\mathbf{X}}^i \right) \frac{1}{\Delta x} }  }{ g_{\mathbf{X}}^{i+1} - g_{\mathbf{X}}^i }.
\end{equation}
Recall that for a given $ \mathbf{x} \in [ x_i, x_{i+1}]$ the conditional distribution function along the second axis is then obtained by
\begin{equation}
 g_{\mathbf{X}=\mathbf{x},\mathbf{Y}}(y) = \frac{g_{\mathbf{X,Y}}}{\int_{ y_{\min}}^{y_{\max}} g_{\mathbf{X,Y}}(\mathbf{x},y)~\mathrm{d}y}
 =  \frac{g_{\mathbf{X,Y}}( \mathbf{x} ,y) }{g_{\mathbf{X}}( \mathbf{x}) },
\end{equation}
where $g_{\mathbf{X}}( \mathbf{x})$ is merely a normalization.
The cumulative conditional distribution function along the second axis reads
\begin{multline}
G_{\mathbf{X}=\mathbf{x},\mathbf{Y}}(y) =   \Bigg\{
\begin{bmatrix}
 1- \frac{ \mathbf{x} - x_{i,j}}{\Delta x}  & \frac{ \mathbf{x} - x_{i,j}}{\Delta x}
\end{bmatrix}
\begin{bmatrix}
 g^{i,j}  & g^{i,j+1}\\
 g^{i+1,j} & g^{i+1,j+1}\\
\end{bmatrix}

\begin{bmatrix}
 \frac{ y - y_{i,j}}{\Delta y}  - \frac{1}{2} \left( \frac{ y - y_{i,j}}{\Delta y} \right)^2 & \frac{1}{2} \left( \frac{ y - y_{i,j}}{\Delta y} \right)^2

\end{bmatrix} \Delta y\\
+
\begin{bmatrix}
 1- \frac{ \mathbf{x} - x_{i,j}}{\Delta x}  & \frac{ \mathbf{x} - x_{i,j}}{\Delta x}
\end{bmatrix}
\begin{bmatrix}
  \sum_{1 \leq k < j}  (g^{i,k}+ g^{i,k+1} ) \frac{\Delta y}{2}\\
  \sum_{1 \leq k < j}  (g^{i+1,k}+ g^{i+1,k+1}) \frac{\Delta y}{2}
\end{bmatrix} \Bigg\} \frac{1}{g_{\mathbf{X}}( \mathbf{x})}

\end{multline}
Once again, given the second uniform random number $ \mathbf{u}_y \in [0,1]$ and the first sample $\mathbf{x}$, we search for the sample $\mathbf{y}$ in the second direction such that
\begin{equation}
G_{\mathbf{X}=\mathbf{x},\mathbf{Y}}(\mathbf{y}) = \mathbf{u}_y.
\end{equation}
Here the cell index $j$ of $\mathbf{y}$ is found by determining the largest $j$ such that
\begin{equation}
\sum_{1 \leq k < j} \left\{
\begin{bmatrix}
 1- \frac{ \mathbf{x} - x_{i,j}}{\Delta x}  & \frac{ \mathbf{x} - x_{i,j}}{\Delta x}
\end{bmatrix}
\begin{bmatrix}
    (g^{i,k}+ g^{i,k+1} ) \frac{\Delta y}{2}\\
    (g^{i+1,k}+ g^{i+1,k+1}) \frac{\Delta y}{2}
\end{bmatrix} \right\}
  <  g_{\mathbf{X}}( \mathbf{x}) \mathbf{u}_y,
\end{equation}
where $i$ still denotes the cell index in the first dimension, $\mathbf{x} \in [x_i, x_{i+1}]$.
Modifying eqn.~\eqref{eqn:bilinear:cellinversionx} by scaling the density with the normalization and interpolation in the first dimension the $y$ coordinate reads
\begin{equation}
\begin{split}
 \mathbf{y} &= y_j + \Delta y
 \frac{ - \gamma_0 + \sqrt{ ( \gamma_0 )^2 + 2 (\gamma_1 -\gamma_0 )
\left(g_{\mathbf{X}}( \mathbf{x}) \mathbf{u}_x - \delta \right) \frac{1}{\Delta y} }  }{ \gamma_1 -\gamma_0 }\\
\text{with}\\
\gamma_l &= \begin{bmatrix}
        1- \frac{ \mathbf{x} - x_{i,j}}{\Delta x}  & \frac{ \mathbf{x} - x_{i,j}}{\Delta x}
      \end{bmatrix}
    \begin{bmatrix}
 g^{i,j+l} \\
 g^{i+1,j+l}
\end{bmatrix},
    \quad l=0,1      \\
\delta&= \sum_{1 \leq k < j} \left\{
\begin{bmatrix}
 1- \frac{ \mathbf{x} - x_{i,j}}{\Delta x}  & \frac{ \mathbf{x} - x_{i,j}}{\Delta x}
\end{bmatrix}
\begin{bmatrix}
    (g^{i,k}+ g^{i,k+1} ) \frac{\Delta y}{2}\\
    (g^{i+1,k}+ g^{i+1,k+1}) \frac{\Delta y}{2}
\end{bmatrix} \right\}.
\end{split}
\end{equation}
\subsection{(Bi)Linear Orthogonal Series Density Estimation}
In order to reconstruct a density from given samples we use the counterpart to the bilinear inverse transform sampling, orthogonal series density estimation with
linear splines. Our linear spline basis $(N^1_j)_{j=1,\dots}$ on a one dimensional grid $(x_j)_{j=1,\dots}$ is defined as
\begin{equation}
\label{linearsplinebasis}
N_i^1(x) = N^1\left(  \frac{x-x_i}{\Delta x} \right), \quad
N^1(x) = \begin{cases}
       1-|x| & \text{ for } |x| \leq 1,\\
           0 & \text{ otherwise}.
          \end{cases}
\end{equation}
Note that the mass matrix required for the $L^2$ projection onto our linear spline basis is sparse has the following coefficients
\begin{equation}
 M^x_{m,n}=  \int N_{m}^1(x) N_{n}^1(x) \dx = \Delta x \begin{cases}
                                    \frac{2}{3}       &\text{ for } m=n,\\
                                    \frac{1}{3}      &\text{ for } |m-n|=1,\\
                                      0    &\text{ otherwise }.
                                  \end{cases}
\end{equation}
For periodic domains the mass matrix is circulant and can be easily applied as an inverse by the use of the fast Fourier transform~\cite{gray2006toeplitz}.
Here, of course, we consider a two dimensional basis such that the bilinear interpolation can be written with the help of tensor product splines as
\begin{equation}
 g(x,y) = \sum_{i,j} g^{i,j} N_{i,j}(x,y), \quad N_{i,j}(x,y):= N^1\left(  \frac{x-x_i}{\Delta x} \right)N^1\left(  \frac{y-y_j}{\Delta y} \right).
\end{equation}
The corresponding mass matrix for the two dimensional space is obtained as a tensor product $M^{x,y} = M^x \otimes M^y$.
This also means that the $L^2$ projection of a density sampled by bilinear inverse transform sampling onto this bilinear spline space is exact:
\begin{equation}
 \kvec{ g_{i,j} } = \left( M^{x,y} \right)^{-1} \kvec{ \iint g(x,y) N_{i,j}(x,y)  \mathrm{d}x \mathrm{d}{y}          }
\end{equation}
The coefficients $g_{i,j}$ can be approximated as
\begin{equation}
 \kvec{ g_{i,j} } \approx \left( M^{x,y} \right)^{-1} \kvec{ \ \frac{1}{N_s} \sum_{n=1}^{N_s} w_n ~ N_{i,j}(\mathbf{x}_n,\mathbf{y}_n) }
\end{equation}
which corresponds to the cloud in cell scheme combined with the mass matrix. Since we are interested in the sampling density the weights are constant $w_n=1$, but for any other
function they read $w_n= \frac{f(\mathbf{x}_n,\mathbf{y}_n)}{g(\mathbf{x}_n,\mathbf{y}_n)}$.\\
Since we combine this OSDE with spectral densities there is an additional approximation error that we have to review.
Recall that the $m^{\mathrm{th}}$ order B-Spline $N^1(x)$ on a grid of size $h$ is obtained by convolution as
\begin{equation}
\begin{split}
N^m(x) &= \underbrace{ N^0 * \dots * N^0}_{m+1 \text{ times }}(x) = (N^0)^{*(m+1)}(x)=\int_{-\infty}^{\infty} N^0(x-y) N^{m-1}(y) ~\mathrm{d}y\\
\text{ with }
N^0 (x) &= \begin{cases}
        \frac{1}{h} &\text{if } x \in \left(-\frac{h}{2}, \frac{h}{2} \right)\\

        0 & \text{else}.
       \end{cases}\\
\end{split}
\end{equation}
The Fourier transform of one basis function $N^m$ reads
\begin{equation}
 \int_{x_{\min}}^{x_{\max}} N^m(x) \euler^{-i k x} ~\mathrm{d}x  = \left[ \mathrm{sinc} \left( \frac{kh}{2} \right) \right]^{m+1},
\end{equation}
such that we can conclude that the relative error on the $k^{\mathrm{th}}$ Fourier mode is
\begin{equation}
 \left| 1- \left[ \mathrm{sinc} \left( \frac{kh}{2} \right) \right]^{m+1} \right|
\end{equation}
For our linear splines $m=1$, this can be quite large on the highest mode still represented on the grid $ 1-\mathrm{sinc}(1)^2 \approx  0.2919$, such that we have to decrease the grid size $h$.
Given the Fourier coefficients for the spectral solver this is easily achieved with increasing the resolution by a factor $N_{\mathrm{pad}}$ by zero padding in the inverse FFT.
When applying padding we chose $N_{\mathrm{pad}}=32$, such that the relative error on highest mode reduces to $ 1-\mathrm{sinc}\left( \frac{1}{N_{\mathrm{pad}}} \right)^2 \approx  3.2548e-04$.
\FloatBarrier
\section{Numerical results}
We consider the two dimensional Vlasov--Poisson system with the initial condition
\begin{equation}
\begin{split}
 f(x,v,t=0) &= \frac{1- \epsilon \cos(k  x)}{\sqrt{2\pi}} \left[  (1-n_b)  \euler^{-\frac{v^2}{2}} + \frac{ n_b}{\sigma_b}  \euler^{-\frac{(v-v_b)^2}{2 \sigma_b^2}}      \right] \\
 & \text{ for } x \in [0,L], ~v \in [v_{\min},v_{\max} ], ~ L=\frac{2\pi}{k}\\
\end{split}
\end{equation}
and the following parameters for two nonlinear test cases:
\begin{center}
\begin{tabular}{c c}

 Landau damping & $ \epsilon= 0.5 , ~ k=0.5, ~[\sigma_b =1 , ~n_b=0]$\\
 Bump-on-tail instability & $ \epsilon= 10^{-3} , ~ k=0.3, ~\sigma_b =0.3 , ~n_b=0.1$\\
\end{tabular}
\end{center}
The number of particles is in the PIC simulation is denoted by $N_p$ and the number of cells by $N_f$. For the Poisson solver finite elements based on cubic B-splines are used.
Random numbers and the quasi random Sobol sequence are provided by MATLAB~\cite{bratley1988algorithm}.
The spectral solver uses the same number of grid points in spatial and velocity space,~$N_x=N_v$.

\FloatBarrier
\subsection{Phase space conservation}
We have already seen in fig.\ref{fig:pic:star_discrepancy} that phase space conserving symplectic Euler as opposed to the explicit Euler preserves the uniformity of a quasi-random sequence implying better convergence.
In general for integrators, which do not preserve phase space volume but are dissipative such as asymptotically preserving schemes
like~\cite{filbet2016asymptotically,filbet2017asymptotically} the likelihoods have to be propagated accordingly. For the explicit Euler, there are two options.
As we are following the characteristics we leave $f_k$ constant, but rescale the sampling likelihood $g_k$ with the according Jacobi determinant of the flow and call this \textit{euler}.
For \textit{euler2} we ignore the characteristics and also rescale $f_k$ with the Jacobian. Since the likelihoods actually change now it makes sense
to look at otherwise conserved quantities such as the total mass
\begin{equation}
 \iint f(x,v,t)  \dxdv = \mathbb{E}\left[ \frac{f(\mathbf{X}(t),\mathbf{V}(t),t)}{ g(\mathbf{X}(t),\mathbf{V}(t),t)}  \right] =
 \frac{1}{N_p} \sum_{k=1}^{N_p} \frac{ f^t_k}{ g^t_k}
\end{equation}
and the discrete variant of the differential entropy

\begin{multline}
\label{discreteentropy}
  \iint f(x,v,t) \ln(f(x,v,t))  \dxdv\\
  = \mathbb{E}\left[ \frac{f(\mathbf{X}(t),\mathbf{V}(t),t) \ln \left( f(\mathbf{X}(t),\mathbf{V}(t),t)  \right)}{ g(\mathbf{X}(t),\mathbf{V}(t),t)}  \right]
  \approx \frac{1}{N_p} \sum_{k=1}^{N_p} \frac{ f^t_k \ln(f^t_k)}{ g^t_k}.
  \end{multline}

Note, that there are different ways to estimate the entropy from a sample,~\cite{beirlant1997nonparametric} gives an overview, while in ~\cite{ebrahimi1994two}
mesh based examples ready for implementation can be found. By propagating the sampling weight correctly it is possible to observe changes in the entropy caused by the dissipative integrator.
The results can be seen for strong Landau damping in fig.~\ref{fig:phasespace_conservation:landau_strong} and the Bump-on-tail instability in fig.~\ref{fig:phasespace_conservation:bumpontail}.
For the electrostatic field energy a reference solution was calculated using a pseudo spectral solver.

\begin{figure}
\centering
\begin{subfigure}{0.8\textwidth}
 \includegraphics[width=\textwidth]{plots/phasespace_conservation/bumpontail_fieldenergy.png}
 \caption{electrostatic energy}
\end{subfigure}\\
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{plots/phasespace_conservation/bumpontail_entropy.png}
\caption{discrete entropy}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
 \includegraphics[width=\textwidth]{plots/phasespace_conservation/bumpontail_energy_error.png}
 \caption{relative energy error}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
 \includegraphics[width=\textwidth]{plots/phasespace_conservation/bumpontail_mass_error.png}
 \caption{relative mass error}
\end{subfigure}\\

\caption{PIC simulation using cubic B-Splines $(N_x=32)$ and $\Delta t=0.01$ of a Bump-on-tail instability with $N_p=10^6$
         quasi-randomly distributed particles drawn by inverse transform importance sampling.
         The explicit Euler (\textit{euler}) is outperformed by the sympletic (\textit{seuler)}).
          }
\label{fig:phasespace_conservation:bumpontail}
\end{figure}

\begin{figure}{H}
\centering
\begin{subfigure}{0.8\textwidth}
 \includegraphics[width=\textwidth]{plots/phasespace_conservation/landau_strong_fieldenergy.png}
 \caption{electrostatic energy}
\end{subfigure}\\
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{plots/phasespace_conservation/landau_strong_entropy.png}
\caption{discrete entropy}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
 \includegraphics[width=\textwidth]{plots/phasespace_conservation/landau_strong_energy_error.png}
 \caption{relative energy error}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
 \includegraphics[width=\textwidth]{plots/phasespace_conservation/landau_strong_mass_error.png}
 \caption{relative mass error}
\end{subfigure}\\

\caption{PIC simulation of nonlinear strong Landau damping with $N_p=10^6$ quasi-randomly distributed particles drawn by inverse transform importance sampling.
         Sympletic (\textit{seuler)})
          Cubic B-Spline $(N_x=32)$ }
\label{fig:phasespace_conservation:landau_strong}
\end{figure}
The choices made for \textit{euler} and \textit{euler2} are both inconsistent in some way because the method itself is just not suited for this purpose.
For an asymptotically preserving scheme, where we do not follow the characteristics, it would make sense though to rescale $f_k$ and $g_k$ accordingly.

\FloatBarrier
\subsection{Convergence of OSDE}
The first experiment should test whether the sampler has been correctly implemented and whether we obtain the expected convergence rates. For this
the Vlasov--Poisson system~\eqref{vlasov1}-\eqref{poisson1} is simulated with the spectral solver for the Bump-on-tail instability and Landau damping in $t \in [0,t_{\max}]$.
The density $f(x,v,t_{\max})$ is obtained on the $N_x \times N_v$ phase-space grid without additional padding. In the following bilinear interpolation
is used for further representation of $f$. Negative values remain and are not truncated.
After obtaining the corresponding sampling density according to eqn.~\eqref{normalizef} samples are drawn by bilinear inverse transform sampling. To verify
that the particles are sampled correctly the density is estimated again by bilinear OSDE.
The convergence rates in fig.~\ref{fig:bumpontail_osde} and fig.~\ref{fig:landau_osde} are as expected such that we can proceed with the PIC simulation.

\begin{figure}{h}
\centering
\begin{subfigure}{0.45\textwidth}
 \includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_unirnd_mc.png}
 \caption{uniform (pseudo)random numbers}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_unirnd_qmc.png}
\caption{uniform Sobol numbers}
\end{subfigure}\\
\begin{subfigure}{0.45\textwidth}
 \includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_samples_mc.png}
 \caption{Random samples}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_samples_qmc.png}
\caption{Quasi-Random samples}
\end{subfigure}\\
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_g.png}
\caption{sampling density $g$ at $t=30$}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_osde_l2_error.png}
\caption{Relative $L^2$ sampling error}
\end{subfigure}
\caption{Sampling $N_p=4000$ random and quasi-random particles with bilinear inverse transform sampling from the sampling density $g$ obtained at $t_{\max}=35$ from spectral simulation
          of the Bump-on-tail instability $N_x=N_v=64,~ \Delta t=0.01$.
          Fig. (f) shows the relative $L^2$ error on the density $f$ obtained by OSDE from the samples at increasing number of particles.}
\label{fig:bumpontail_osde}
\end{figure}

\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
 \includegraphics[width=\textwidth]{plots/it_sampling/landau_512_samples_mc.png}
  \caption{Random samples}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/landau_512_samples_qmc.png}
\caption{Quasi-Random samples}
\end{subfigure}\\
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/landau_512_g.png}
\caption{sampling density $g$ at $t=30$}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/landau_512_osde_l2_error.png}
\caption{Relative $L^2$ sampling error}
\end{subfigure}
\caption{The fine structures of nonlinear Landau damping are resolved at $t_{\max}=30$ by a spectral solver with $N_x=N_v=512,~ \Delta t=0.01$ in (c). The difference between
         $N_p=4000$ random (a) and quasi-random (b) samples obtained from (c) are clearly visible and also seen in (d) in the different convergence rates
         of the corresponding reconstructing OSDE of $f$.}
\label{fig:landau_osde}
\end{figure}

\FloatBarrier
\subsection{Density Estimation by Bilinear Interpolation}
Using a Monte Carlo based density estimator for recovering the density $f$ from a distribution of markers is only using the information $\frac{f}{g}$ and neglecting the fact that each marker
already transports the value of the density $f$ or $g$ respectively. Using this additional information leads
to an interpolation problem. We sampled from a bilinear interpolant such that it is reasonable to test the least square fit of the bilinear interpolation coefficients under different marker distributions.
If there are less more grid points than markers the problem is well-posed, but otherwise one has to add a
regularization to the least square problem. The easiest choice was $L^2$ regularization known as ridge regression. Figure~\ref{fig:bumpontail_interp} shows as expected that interpolation delivers better results than Monte Carlo OSDE. For a high precision reconstruction interpolation relies on the uniformity of the interpolation points, such that uniform sampling outperforms the importance sampling by magnitudes.
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
 \includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_nl_128_f_interp_conv.png}
 \caption{Convergence of reconstruction}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_nl_128_f_interp_qmc.png}
 \caption{Bilinear reconstruction of $f$}
\end{subfigure}\\
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_nl_128_samples_qmc_interp.png}
\caption{Quasi-Random samples}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/it_sampling/bumpontail_nl_128_f_interp_qmc_l1error.png}
\caption{absolute error on  of $f$}
\end{subfigure}

\caption{The density of the Bump-on-tail instability calculated by a spectral solver $N_x=N_v=128,~ \Delta t=0.01$ at $t=100$ is sampled with $N_p=10^4$ quasi-random particles in (c). Figure (b) shows the reconstruction
from these markers by bilinear interpolation with ridge regression. Due to the importance sampling
the error (d) is dominant in the regions of low but nonzero density. In (a) the error of the reconstruction by interpolation is shown for a different number of samples, where the break-even point lies at the number of degrees of freedom $N_x N_v$. When there are fewer markers than grid points and the problem is ill-posed importance sampling performs best but it does not converge for the well-posed problem. Uniform sampling is only better in the well-posed regime, where it converges to machine precision. It becomes also clear that the Quasi-Random points are more uniform than the random ones.}
\label{fig:bumpontail_interp}
\end{figure}{h}

\FloatBarrier
\subsection{Eulerian and Lagrangian Coupling}
Since we are confident that particles can be sampled correctly, we proceed with the coupling between the Eulerian and Lagrangian solver. Coupling Eulerian and PIC codes is a rapidly developing field~\cite{dominski2018tight,choi2018coupling}.
First the Bump-on-tail instability is considered, which starts with a very small amplitude posing no problem for the spectral
solver $(N_x=N_v=32,~v_{\max}=-v_{\min}=10,~\Delta t=0.1)$, which follows the reference solution $(N_x=N_v=8 \times 32= 512)$ properly in the linear phase~\ref{fig:coupling:bumpontail:fieldenergy}
but later suffers from oscillations due to filamentations and aliasing. Therefore, we switch to PIC ($N_p=10^6~,N_f=16$) at $t_0=35$.
Since the resolution is quite low and we want to suppress aliasing due to the low order interpolation the spectral density is zero padded with a factor of $N_{\mathrm{pad}}=32$ yielding
a $1024 \times 1024$ phase space grid for the sampler. Given the curse of dimensionality high order splines would be a better solution.
Nevertheless the PIC code follows the reference visibly better in fig.~\ref{fig:coupling:bumpontail:kineticenergy} and fig.~\ref{fig:coupling:bumpontail:kineticenergy}.
The third order symplectic Runge Kutta time discretization is the same for PIC and the spectral solver, but
PIC is derived from a Lagrangian formulation such that the better energy conservation in fig.~\ref{fig:coupling:bumpontail:energy_error} is no surprise.\\
Although fig.~\ref{fig:bumpontail_osde} and fig.~\ref{fig:landau_osde} confirm the superior convergence rate of the QMC numbers for the initial sampling, the question remains whether this stays true
over the nonlinear phase. This question was already addressed before~\cite{ameres2016adv,ameres2018} and can also be answered positively here by fig.~\ref{fig:coupling:convgc:bumpontail}.
Note that the symplectic Runge Kutta scheme preserves phase space volume, hence the Jacobian of the discrete flux is exactly one which preserves also the likelihood of each marker.
This means that the discrete flux induces a measure-preserving map such that the Hausdorff measure is preserved. Therefore, the uniformity of the low discrepancy sequence is preserved such that
the higher order convergence rates for QMC keep their validity~\cite{basu2016transformations,aistleitner2013low,de2018quasi}.

The only confusing issue that can emerge is, that the total variation (QMC) as well as the variance (MC)
of the \textit{entire} map from the initial condition to a certain time $t$ increases with the nonlinearities of the transport. This, however, does not change the convergence rates at a certain time $t$.\\
The same procedure with $(N_x=N_v=64,~v_{\max}=-v_{\min}=6.5,~\Delta t=0.05,~t_0=30,~N_p=10^6,~N_f=16, N_{\mathrm{pad}}=32$ is repeated for nonlinear Landau damping.
Here, because of the many perturbations at small amplitudes the difference between the spectral solver and PIC remains small, see fig.~\ref{fig:coupling:landau}.
\begin{figure}
\centering
\begin{subfigure}{0.7\textwidth}
\includegraphics[width=\textwidth]{plots/spectral_pic_coupling/bumpontailfieldenergy.png}
\caption{electrostatic energy}
\label{fig:coupling:bumpontail:fieldenergy}
\end{subfigure}\\
\begin{subfigure}{0.45\textwidth}
 \includegraphics[width=\textwidth]{plots/spectral_pic_coupling/bumpontailenergy_error.png}
 \caption{relative energy error}
 \label{fig:coupling:bumpontail:energy_error}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/spectral_pic_coupling/bumpontailkineticenergy.png}
\caption{kinetic energy}
\label{fig:coupling:bumpontail:kineticenergy}
\end{subfigure}
\caption{Transition from a spectral solver to PIC for the Bump-on-tail instability. The PIC code is initialized at $t_0=35$ by QMC inverse transform sampling.
         Although it appears from (a) and (b) as if the spectral solver is noisier than the PIC code this effect caused by filamentation is suppressed by better resolution in the reference.
         Depending on the implementation, the spectral solver, in general, outperforms the PIC code in two dimensions $(32^2<512^2<10^6)$.}
\label{fig:coupling:bumpontail}
\end{figure}

\begin{figure}
\centering
 \includegraphics[width=0.7\textwidth]{plots/spectral_pic_coupling/bumpontail_fieldenergy_convgc.png} -
\caption{$L^2$ error of the electrostatic energy for quasi random (QMC) and random (MC) numbers in the spectral-PIC coupling for bump-on-tail instability from fig.~\ref{fig:coupling:bumpontail}.
         In order to dominantly include nonlinear effects the error is only taken for $t \in (39.9, 49.9)$. Here QMC performs better even with the included nonlinear effects and the excluded initial sampling.}
 \label{fig:coupling:convgc:bumpontail}
\end{figure}

\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
 \includegraphics[width=\textwidth]{plots/spectral_pic_coupling/landau_fieldenergy.png}
\caption{electrostatic energy}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{plots/spectral_pic_coupling/landau_energy_error.png}
 \caption{relative energy error}
 \end{subfigure}

\caption{Transition from a spectral solver to PIC for nonlinear Landau damping. The transition at $t_0=30$ is precisely chosen in moment of small amplitudes in order to raise the difficulty for the PIC.}
\label{fig:coupling:landau}
\end{figure}

\FloatBarrier
\section{Conclusion and Outlook}
\begin{figure}[h]
\centering
\includegraphics[width=0.75\textwidth]{plots/it_sampling/eulerianlagrangianboundary.png}
 \caption{The image shows the upper half of the phase space of a Vlasov--Poisson simulation with a transition from a Eulerian to a Lagrangian region. This means the entire phase-space is flowing to the right. The white dashed line marks the boundary/interface position. During the split-step $\dot{X}=V$ the red triangle encloses the area of phase space entering the Lagrangian domain at each time step of size $\Delta t$. In this case, the Eulerian solver can also precalculate the volume in the red triangle, such that it can be sampled directly via inverse transform sampling at each time step, which results in a smooth density. If the boundary condition is only given at a fixed point $\bar{x}$ with density
 $f_{\mathbf{X}=\bar{x}}(v)$ it is not sufficient to draw the particles velocity as $\mathbf{V} \sim f_{\mathbf{X}=\bar{x}} $ and set $\mathbf{X}=\bar{x}$. To account for a volume one has to draw
 $\mathbf{X} \sim \mathcal{U}(\bar{x}, \bar{x} + \Delta t \mathbf{X})  $ uniformly. For a strong electric field, the splitting might not be appropriate such that the phase space integral becomes more complicated but the principle stays the same.}
\label{fig:eulerianlagrangianboundary}
\end{figure}
It was shown how to use Quasi Monte Carlo numbers to sample from arbitrary phase-space densities. Not only does this improve convergence but opens up the possibility to couple established Eulerian codes to the big PIC codes for kinetic and gyrokinetic systems.
It should also be noted that
Considering that Vlasov problems are high dimensional one has to realize that computing high dimensional marginals can be expensive.
Here~\cite{olver2013fast} proposed to use a low-rank spectral representation, which is a major underlying idea of Approxfun.jl~\cite{olver2014practical}.
For the Vlasov--Poisson system there exist solvers in tensor train format~\cite{kormann2015solving} that already provide a low-rank approximation of $f$
which possibly can be exploited for efficient sampling in the future.
Another aspect not followed here is the transition from a Lagrangian to a Eulerian solver. OSDE is for sure a possibility but this misses the point of Lagrangian solvers completely such that an interpolation approach as shown in fig.~\ref{fig:bumpontail_interp} is recommended.
Each particle transports the value of the density along the characteristics such that the key for such a transition lies in a suitable interpolation, which has already been extensively discussed for Semi-Lagrangian solvers~\cite{sonnendrucker1999semi,bermejo1991analysis}.
In two dimensions phase-space conservation is the same as symplecticity, but in higher dimensions, symplecticity is something mildly stronger. Hence, it remains the question, what impact symplecticity has on the discrepancy.\\
Another important application of the sampling techniques presented here is the implementation of boundary conditions for the Vlasov equation in PIC. The standard approach is to draw particles at an interface position according to a velocity distribution at each time step.  PIC relies on phase space conservation such that any insertion of markers should actually be an insertion of a volume of phase space. This means that the time step cannot be assumed to be infinitesimally small, but also the boundary condition has to be integrated exactly over time. In most cases, this can be done analytically and otherwise, a numerical pre-calculation is sufficient. As explained in fig.~\ref{fig:eulerianlagrangianboundary} this requires sampling nontrivial phase space volumes, where this article opens new possibilities.

\section{Acknowledgement}
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

\FloatBarrier

\printbibliography

\end{document}

config: arxiv


TITLE: Prove a complex number is real
QUESTION [4 upvotes]: Let $z$ be a complex such that $|z-1| =1$, and consider the complex numbers $v$ and $w$ such that: $w = z^2 -z$ and $3\arg(v) = 2\arg(w)$, where arg  is the argument of a complex number. Show that $$ \left({zv\over w}\right)^2$$ is a real number
I first tried to show that the imaginary part of the number is zero, knowing that half of the subtraction of a complex number and its conjugate is equal to its imaginary part. Nonetheless, i could not find a relation between $z$ and $v$, and I ended up in a answer with these two letters.

REPLY [1 votes]: The answer for two non real numbers is possible as $i^2=-1$ but it isnt the case for $i^3=-i$ so we cant prove it for three non real numbers also if for a quadratic if roots are imaginary then they exist as conjugates and their product is a real number. Hope this helps you.Let the three numbers be $1+0i,2+0i,-3+0i$ so we directly get the result . I dont think there exists any rigorous way to prove it.REAL numbers are subset of complex numbers hence I have given you this example.


TITLE: Is there a name for this identity-like object?
QUESTION [3 upvotes]: I was grading for a linear algebra class just now and someone remarked that since $U + U = U$, then $U$ would be an identity for the set of subspaces of some vector space with binary operation subspace addition. Obviously I corrected their mistake, since it isn't true that this fixed $U$ has the property $U + V = V$ for all subspaces $V$ of this vector space (except perhaps in trivial cases), but it made me curious whether or not there is a name for this type of object since they have exhibited a nontrivial instance of the existence of such an element.
I suppose it would only be meaningful in a semigroup or algebraic object with less structure than a semigroup, since the standard proof that the identity element is unique in, say, a group, would probably show that any element like this would have to be the identity.
So, my question is, is there a name for an element of a semigroup (or an algebraic object with less structure if even a semigroup in general is too restrictive), let's call the element $x$ and the operation of the semigroup +, with the property that $x + x = x$.

REPLY [3 votes]: This property is called idempotence, and it can be nontrivial even in monoids; for example, there is a very large class of monoids all of whose elements are idempotent, namely any semilattice. In particular, the collection of subspaces of a vector space is an example. There are two corresponding idempotent monoid operations, namely intersection and sum.
Idempotents behave like identities in the following sense: if $m$ is an idempotent element of a semigroup $S$, you can consider the subsemigroup of elements $a$ satisfying $ma = am = a$. This is also a semigroup, even a monoid, and (by construction) has identity element $m$; in fact by construction it's the largest subsemigroup of $S$ on which $m$ acts as the identity.
When applied to an element $U$ of the lattice of subspaces of a vector space $V$, with monoid operation given by intersection, this construction produces the lattice of subspaces of $U$. With the monoid operation given by sum, this construction produces the lattice of subspaces of $V$ containing $U$, which can naturally be identified with the lattice of subspaces of $V/U$ by the isomorphism theorems.


TITLE: Magnetic induction: apparent contradiction between local and integral prediction
QUESTION [1 upvotes]: I'm trying to understand better magnetic induction, and I came across a situation where I find an apparent contradiction between the local form ($\nabla×E=-\partial B/\partial t $) and the corresponding integral form ($e=-d\phi/dt$, which of course corresponds to the integration of the local form) of induction law.
[EDIT : in order to make the problem easier to visualize in space, I found a simpler situation with a similar problem]
Take these two coaxial coils :

We provide an increasing current $i_1$ to coil 1, and want to examine the induced current $i_2$ in coil 2.
The question is : in which direction does $i_2$ flow?
Integral reasoning
Surface of reference : the cross section of coil 2 (in green).

The magnetic field $B_1in$ in coil 1 is oriented up.
There's also a downward-pointing magnetic field $B_1out$ outside coil 1, in the small annular zone between coil 2 and coil 1.
But globally the integral $\phi$ of the field over S is positive.

So, as $i_1$ increases, $\phi$ is more and more positive.  Therefore $e=-d\phi/dt <0$ : $i_2$ is such that the induced field points downward, in opposition to the increase of $\phi$ ($i_2$ < 0)
Local reasoning
the wire of coil 2 is subject to $B_1out$ : the wire doesn't "know" what happens inside coil 1. If, instead of coil 1, we had an externally-imposed flux pointing downward everywhere, then coil 2 couldn't tell the difference.
As $i_1$ increases, $B_1out$ gets stronger, so $\nabla×E=-\partial B/\partial t $ points up at each point of the wire loop. If the curl points up, it means that the induced electric field is such that $i_2$ is positive (the field induced in coil 2 points up).
This is obviously a contradiction with the integral reasoning. I know that the correct solution is the first one, but I don't see where the mistake is in the local reasoning.
[old description]
The situation is simple: a permanent magnet rotates counter-clockwise (the axis of rotation is z) in a wire loop (the loop is in the y-z plan).

I want to examine the induced current at instant t when the magnet "crosses" the plan of the loop (that is to say : the moment seen on the picture above).
If I use the integral form of Faraday's law, I have an integrated flux that is positive at t-dt,  and negative at t+dt. So, the variation of the flux is negative, and the induced current is positive. The induced magnetic field points in the +x direction (the induced magnetic field opposes the rotation of the permanent magnet, which is expected).

So far, so good.
Now, let's have a look at the magnetic field that is applied to the wire, to use local form of the equation.
If we look from above, the magnetic field at t-dt points slightly to the left. At t+dt, it points slightly to the right (it's crucial to notice that the field lines aren't quite parallel, and that they curve out of the magnet axis faster that the magnet rotates). So, $\partial B/\partial t $ has a positive x component, which means that $\nabla \times E$ points to the left (negative x component)... which is in total contradiction with the prediction of the integral form.

We have the same problem if we look at the horizontal portions of the loop: the magnetic field points in the opposite direction of that of the magnet, due to the fact that the magnetic lines form a closed loop.

So, when the magnet rotates, the magnetic field variation in the horizontal wire points to the right, and the curl of the induced current points to the left (as well as the induced magnetic field).

You can drive the same analysis at any point of the wire : the induced magnetic field points to the left, because the magnetic field outside the magnet actually isn't aligned with that "inside" the magnet.
We see therefore that the integral form predicts something different from the local examination.
I know that the local analysis I made is wrong, because it predicts that the induced field helps the magnet to rotate (so, any slight pulse at the beginning would lead to a never-ending self-acceleration of the magnet, creating mechanical energy from nowhere).
Yet, I don't see what the flaw is in my local analysis. I mean : the wire actually is subject to the local magnetic field, right? And this magnetic field variation points towards the positive x, so...
I don't get what I'm missing.

REPLY [0 votes]: Haven't read all of this, But you confusion seems to stem from the orientation of the surface.
If there is a surface that the orientation is positive on the top.
If there is a magnetic flux through this that the field points upwards
Then da and B are in the same direction, so there is a positive flux. If this flux increases then dphi/dt is positive , so so the line integral of E .da would be negative ( as -positive) , the line integral in this manner is defined anticlockwise. meaning negative work is done going anticlockwise. aka the field then goes clockwise
The differential form. if db/dt is positive then the curl is negative, ( as -db/dt) if the curl is negative then the field goes clockwise relative to my orientation
there is no  contradiction.
stokes theorem is only valid, if the line integral is taken to be the anticlockwise path, Given my orientation is "up", it is defined this way because the curl is positive, if there is positive "circulation" going anticlockwise ( relative to the curl vector bein taken as "positive") thus positive curl dot positive da, gives me a positive number , this clearly means that stokes theorem would give me a positive number going anti clockwise
If da is orientation is downwards, and the curl is up,  curl E dot da would be negative, meaning relative to this new orientation negative work would be done going anticlockwise, which means the field points clockwise.
Clockwise for the upside down orientation is anticlockwise for the positive orientation