12 of 14 added 3 characters in body

tl;dr

Final results first:

(*Function Definition*)
ClearAll[opaFun];
Options[opaFun] = Options[ListPlot];
opaFun[points_, opts : OptionsPattern[]] := Module[{f, steps = 10 },
   f[x_] := Min[Norm /@ Flatten[ImageData@
         ListPlot[points, opts, Axes -> False, PlotStyle -> {Black, Opacity[x]}]], 1]/Sqrt@3;
   Return@NestWhile[Append[#, {#[[-1, 1]] + 1/steps, f[#[[-1, 1]] + 1/steps]}] &, {{0, 1}},
                    UnsameQ @@ ({##}[[-1, All, 2]]) &, 2, steps + 1][[-2, 1]];];
(*Usage*)
opts = {AspectRatio -> 1/2, ImageSize -> 300, Axes -> False};
Timing@ListPlot[#, opts, PlotStyle -> {Black, Opacity[opaFun[#]]}] &@
                  RandomReal[NormalDistribution[], {10000, 2}]

Mathematica graphics

Now the full answer:

A little but time consuming experiment first. You'll not need to do it, as it is only for finding a model.

Let's see how the Min value of the plot intensity varies with Opacity:

theData = RandomReal[NormalDistribution[], {2000, 2}];
f[x_] := f[x] = 
   Min[Norm /@ Flatten[ImageData@Rasterize[
        ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200, 
                   PlotStyle -> {Black, Opacity[x]}, Axes -> False]], 1]];
Plot[f[x] , {x, 0, .4}, PlotRange -> Full]

Mathematica graphics

So, it is an exponential. (Note: in the last edit to this post I got rid of the exponential model by fitting a few points with an Interpolation, and it works pretty nice)

Let's fit it:

data = Table[{i, f[i]}, {i, 0, 1, .1}]

model = a Exp[b x];
fit = FindFit[data, model, {a, b}, x];
modelf = Function[{t}, Evaluate[model /. fit]]

Show[ListPlot@data, Plot[modelf[x], {x, 0, 1}]]

Mathematica graphics

Now you are ready to set the min value of the brightness of the plot to whatever you want:

(The Sqrt@3 is a normalization factor for the intensity of the {1,1,1} RGB pixel.)

Let's use it:

opac = x /. Solve[# == a E^(b x)/Sqrt@3, x] /. fit & /@ {1/2, 1/4, 1/20, 1/200}

ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200, 
   PlotStyle -> {Black, Opacity[#[[1]]]}, Axes -> False] & /@ opac

Mathematica graphics

Edit

Let's pack that all together in a function and plot two very different point sets with the same maximum darkness.

ClearAll[opa];
Options[opa] = Options[ListPlot];
opa[desiredOpacity_, points_, opts : OptionsPattern[]] :=
 Module[{f, a, b, model, fit, modelf, x},
  f[x_] := 
   f[x] = Min[
     Norm /@ Flatten[
       ImageData@
        Rasterize[
         ListPlot[points, Axes -> False, 
          PlotStyle -> {Black, Opacity[x]}]], 1]];

  model = a Exp[b x];
  fit = FindFit[Table[{i, f[i]}, {i, 0, 1, .1}], model, {a, b}, x];
  modelf = Function[{t}, Evaluate[model /. fit]];
  Return[x /. Quiet@Solve[# == modelf[x]/Sqrt@3, x][[1]] /. fit &@
    desiredOpacity];
  ]

theData  = RandomReal[NormalDistribution[], {2000, 2}];
theData1 = RandomReal[NormalDistribution[], {10000, 2}];

opad = opa[.5,  theData,  AspectRatio -> Automatic, ImageSize -> 200];
opad1 = opa[.5, theData1, AspectRatio -> Automatic, ImageSize -> 200];

Grid[{{ListPlot[theData, Axes -> False, PlotStyle -> {Black, Opacity[opad]},
                         AspectRatio -> Automatic, ImageSize -> 200], 
       ListPlot[theData1,Axes -> False, PlotStyle -> {Black, Opacity[opad1]}, 
                         AspectRatio -> Automatic, ImageSize -> 200]}}, 
   Frame -> All]

Mathematica graphics

The same, but darker

Mathematica graphics

Remember that the default plot is:

Mathematica graphics

Edit

Answering @Oleksandr comments, the following does not assume an exponential model:

ClearAll[opa];
Options[opa] = Options[ListPlot];
opa[desiredOpacity_, points_, opts : OptionsPattern[]] := 
 Module[{f, model, x}, 
  f[x_] := f[x] = 
    Min[Norm /@ Flatten[ImageData@ Rasterize[
         ListPlot[points, Axes -> False, PlotStyle -> {Black, Opacity[x]}]], 1]];
  model = Interpolation[Table[{i, f[i]}, {i, 0, 1, .1}]];
  x /. FindRoot[desiredOpacity == model[x]/Sqrt@3, {x, 0, 1}][[1]]]