Disclaimer: this can be done using gnuplot as mentioned in this answer, but you should probably consider another tool for drawing this particular type of chart.
There is at least one way to do this with data preprocessing. The idea is to simulate a glow effect using a Gaussian core to smear data points. Consider the following data contained in a file called data :
1 2 1 2.1 1.1 2.2 2 3 3 4
I intentionally placed the first 3 points close to each other in order to be able to observe the enhanced glow of neighboring points. This data is as follows:

Now we spread the data points using a two-dimensional Gaussian kernel. I wrote the following python code to help with this. The code has a cut-off of four standard deviations ( sx and sy ) around each point. If you want the glow to be a circle, you must choose standard deviations so that the sx / sy ratio is the same as the x / y axis length ratio in gnuplot. Otherwise, the dots will look like ellipses. This is the code:
import numpy as np import sys filename = str(sys.argv[1]) sx = float(sys.argv[2]) sy = float(sys.argv[3]) def f(x,y,x0,y0,sx,sy): return np.exp(-(x-x0)**2/2./sx**2 -(y-y0)**2/2./sy**2) datafile = open(filename, 'r') data = [] for datapoint in datafile: a, b = datapoint.split() data.append([float(a),float(b)]) xmin = data[0][0] xmax = data[0][0] ymin = data[0][1] ymax = data[0][1] for i in range(1, len(data)): if(data[i][0] < xmin): xmin = data[i][0] if(data[i][0] > xmax): xmax = data[i][0] if(data[i][1] < ymin): ymin = data[i][1] if(data[i][1] > ymax): ymax = data[i][1] xmin -= 4.*sx xmax += 4.*sx ymin -= 4.*sy ymax += 4.*sy dx = (xmax - xmin) / 250. dy = (ymax - ymin) / 250. for i in np.arange(xmin,xmax+dx, dx): for j in np.arange(ymin,ymax+dy, dy): s = 0. for k in range(0, len(data)): d2 = (i - data[k][0])**2 + (j - data[k][1])**2 if( d2 < (4.*sx)**2 + (4.*sy)**2): s += f(i,j,data[k][0],data[k][1],sx,sy) print i, j, s
Used as follows:
python script.py data sx sy
where script.py is the name of the file where the code is located, data is the name of the data file, and sx and sy are standard deviations.
Now, back to gnuplot, we define a palette that mimics a luminous pattern. For isolated points, the summed Gaussians give 1 at the point position; for overlapping points, it gives values above 1. You must take this into account when defining a palette. The following is an example:
set cbrange [0:3] unset colorbox set palette defined (0 "black", 0.5 "blue", 0.75 "cyan", 1 "white", 3 "white") plot "< python script.py data 0.05 0.05" w image

You can see that the points are actually ellipses, because the ratio of the axis lengths does not coincide with the ratio of standard deviations in different directions. This can be easily fixed:
plot "< python script.py data 0.05 0.06" w image
