Generate a point cloud distributed according to a Boltzmann-Gibbs distribution with prescribed marginals Let $p$ be a probability density on $\Omega\in\mathcal B(\mathbb R^d)$ for some $d\in\mathbb N$ (I'm primarily interested in $\Omega=[0,1)^d$). We can approximate $p$ by $$A_x(y):=\sum_{i=1}^k\varphi_{x_i}(y):=\sum_{i=1}^kp(x_i)e^{-\frac{\|x_i-y\|^2}{\sigma(x_i)^2}}$$ for $x\in\Omega^k$ for some $k\in\mathbb N$ and $y\in\mathbb R^d$. $\sigma$ must be suitably defined and $x\in\Omega^k$ be chosen such that $$E(x):=\int_\Omega|A_x(y)-p(y)|\:{\rm d}y$$ is minimized.
Now let $\lambda$ denote the Lebesgue measure on $\mathcal B(\mathbb R)$, $$\varrho(x):=e^{-\frac{E(x)}T}\;\;\;\text{for }x\in\Omega^k$$ for some $T>0$ and $$c:=\int\varrho\:{\rm d}\left.\lambda^{\otimes d}\right|_\Omega^{\otimes k}.$$ Let $\mu$ denote the measure with density $\frac\varrho c$ with respect to $\left.\lambda^{\otimes d}\right|_\Omega^{\otimes k}$ and assume $X$ is a $\Omega^k$-valued random variable distributed according to $\mu$.
If $T\to0+$, then $$\operatorname P\left[X\in B\right]\xrightarrow{T\to0+}\left.\lambda^{\otimes d}\right|_\Omega^{\otimes k}(B\cap\{E=0\})\tag1$$ for all $B\in\mathcal B(\Omega)^{\otimes k}$. This should show that $E(X)=0$ (in the limit $T\to0+$) almost surely and hence $$A_X=p\;\;\;\text{almost everywhere}.$$

However, are (in some sense) the components $X_i$ distributed according to $p$?

Remark/Context: I'm trying to understand the assertions which are made in this paper.

EDIT
Sextus Empiricus mentioned in a comment below his answer that "On average, the density of points will equal $p$". But what does that mean exactly?
Does it mean that the $X_i$ have point density $\pi:=p\left.\lambda^{\otimes d}\right|_\Omega$ (measure with density $p$ with respect to $\left.\lambda^{\otimes d}\right|_\Omega$?
The "intensity" of the point process $(X_1,\ldots,X_k)$ is defined as $$\alpha(B):=\operatorname E\left[\kappa(B)\right]=\sum_{i=1}^k\operatorname P\left[X_i\in B\right]\;\;\;\text{for }B\in\mathcal B(\Omega),$$ where $\kappa$ denotes the random measure $$\kappa(B):=\sum_{i=1}^k\delta_{X_i}.$$
So, the "average number of points in $B\in\mathcal B(\Omega)$" should be $\mu(B)$. Maybe the claim is that $\mu$ has density $p$ with respect to $\left.\lambda^{\otimes d}\right|_\Omega$. But how can we show this?


However, are (in some sense) the components $X_i$ distributed according to $p$?

The image 2 in the paper demonstrates this

The kernels will be taller in places with higher value $p$ but also less wide, and each kernel contains the same mass. To get a higher values of the approximating function in regions where $p(x)$ is high, the points/kernels $x_i$ will need to be closer to each other in places with higher $p$. Note in the image that on the right side of the trapezium, where the green line is higher, the points/kernels are closer to each other.

This shows that $E(X)=0$ (in the limit $T\to0+$)

The case $E(X) = 0$ is not in the sample space. In the limit of the temperature going to zero, it is the lowest energy configuration of the points that will be the most likely sample of points $x$. For the case of $p$ being homogeneously distributed in a two dimensional plane, the lowest energy occurs with a hexagonal pattern of the points.
Example 1
Below is an example toy model that approaches a simplified version of the article that you referred to.


*

*Initiation: I sample 200 points $x_i$ uniformly in a 1 by 1 square $x_i  \in [0,1]^2$

*Function to compute probability: I use a function to compute the 'energy' that computes the difference between a Gaussian mixture model that places a kernel at each point $x_i$ and a function $p$ that uses the sum of coordinates (ie $p=x+y$ such that it is higher in the upper right and lower in the lower left). I integrate the squared difference by using a grid of 100 by 100.

*Algorithm for sampling: I select 5 random points, change their coordinates with a small Gaussian distributed deviation, and compute the 'energy' for this new setup. Then I roll a 'dice' that lands with probability $e^{-E_{old}/T}$ on the old state and with probability $e^{-E_{new}/T}$ on the new state. If the dice lands on the new state, then I change the current state, otherwise I try sampling again.

In the image you see the samples being produced by this procedure. Initially the energy is high (because I started with a uniform distributed sampling) but after some time the points hover around an energy level of 0.2.
Density What you can see in the image is that the points in the upper right corner are closer to each other than the points in the lower right corner. But this is blue noise. If the temperature is very low then some points will be fixed to a particular place. If you look at very small scales, then there might be spaces that are always empty and spaces that are always filled, and not relating to the value $p$. However, when you look at a larger scale, and pick randomly a region, then the number of particles will approximate the density of the function $p$ inside that box.
E = 0 You do not get a state where $E=0$ when $T \to 0$. The only states that are possible are samples of $n$ points in that square, for which the energy will be non-zero. What happens with smaller $T$ is that the states with lower energy (but lower entropy as well) become more likely. Eventually, the lowest energy state (which might possibly be only a single configuration of the points) will be the most likely, and approach probability 1, if the temperature approaches zero.
Additional note to example 1
A confusion might arrise from two spaces that are present. This is shown in the image below.
On the left is the $n$-dimensional space in which a point cloud of $k$ points is drawn (the particular example has $n=2$). Each particular configuration of $k$ points represents itself a single point in a $n \times k$-dimensional space (depicted on the right of the image).
In the first space (on the left) the density of the points approaches $p(x)$ (more or less with some slight details, as the blue noise is not the same as white noise and for small $T$ you get a configuration that favours certain places, see for instance the hexagonal pattern for approximating a uniform $p(x)$).
In the second space (on the right) the density of the points is given by $\varrho(\vec{x}) \propto e^{-E(\vec{x})/T}$. Unless there is a cloud point $\vec{x}$ that can perfectly approximate the function $p(x)$ there are no points in this space for which $E(\vec{x}) = 0$.

The marginal distribution in the right image, the distribution of a coordinate $(x_i)$ for a particular point $i$, can be seen to be related to $p(x)$, because it relates to the points in the left image. The point density in the cloud relates to $p(x)$ and all points are indistinguishable and have equal probability to be in a particular space, so the distribution of a single point does also relate to $p(x)$.
Example 2
Let's consider a simpler case with 2 kernels/points/components giving a function $f(x)$ (or $f(x;x_1,x_2)$ when we consider it parameterized by the two points) to represent a density function $p(x)$ from zero to one. We can represent the 2 points together as a vector $(x_1,x_2) \in [0,1]^2$ and compute the energy and probability for the points. Below is an example of $f(x)$ when it is the sum of two Gaussian components in the points $x_1 = 0.25$ and $x_2 = 0.75$

We can vary the points $x_1$ and $x_2$ and for each case compute the difference $$SS = \int_0^1 (f(x; x_1, x_2) - p(x))^2 dx$$
The mean difference between the Gaussian mixture and the distribution $p(x)$ that we estimate.
Based on $SS$ we can compute the density $$f(x_1,x_2) \propto e^{-SS/T}$$

*

*The image below shows these results. For all possible samples $x_1,x_2$ the energy is non-zero.

*When we decrease the temperature, then the low values $SS$ become more probable.


The plots where created in R with the function
SS = function(x1,x2) {
  sig = 2/sqrt(2*pi)*0.25 ## if p(x) is not uniform, then this sigma might change as function of x2 and x2
  y = 0.5*(dnorm(xi,x1,sig)+dnorm(xi,x2,sig)) # the function f(x; x1,x2)
  mean((y-1)^2) # the mean of the square difference (f(x; x1,x2)-p(x))^2
}

