How can I generate a legend for a BubbleChart?

Since BubbleChart does the scaling of the bubbles the trick is to let it draw the circles of the legend too. I do that by adding the legend circles to the plot data. The rest of the graphics can be added using Show.

lc = Append[Reverse@CityData[#, "Coordinates"], 
            CityData[#, "Population"]
     ] & /@ CityData[{Large, "USA"}];

steps = 8;
legendCircles = Table[{-120 + 18/steps i^1.25, 22, 1000000 i}, {i, steps}];
text = Table[Text[i, {-120 + 18/steps i^1.25, 19}], {i, steps}];

Show[
  Graphics[{FaceForm[LightGray], EdgeForm[Black], 
            CountryData["USA", "Polygon"], text,
            Text[Style["Population (mln)", FontFamily -> "Arial", Bold, 16]]}, 
            PlotRange -> {{-126`, -66}, {17, 51}}, AspectRatio -> Automatic
  ], 
  BubbleChart[Join[lc, legendCircles], BubbleScale -> "Area", 
              AspectRatio -> Automatic, BubbleSizes -> {0.01, 0.1}
  ]
]

or with labels shifted to appear in the bubbles:

---

A note on combining BubbleChart graphics with other graphics

You might have noticed the above use of AspectRatio and PlotRange. They were pretty much essential. Since BubbleChart uses the AspectRatio of its plot to determine the eccentricity of the ovals it draws (yes, it draws ovals: check its FullForm), it is important to watch the aspect ratio of a combined plot carefully.

For instance, let's see what happens if we start from scratch and want to have a cartographic projection different from the simple "Equirectangular" lat,long coordinate system used above.

(* conversion of city locations to Mercator projection *)
lc = Append[
     GeoGridPosition[
       GeoPosition[
        Append[CityData[#, "Coordinates"], 
          CityData[#, "Elevation"]] /. Missing["NotAvailable"] -> 0], 
       "Mercator"][[1, ;; 2]]
     , CityData[#, "Population"]
     ] & /@ CityData[{Large, "USA"}];

(* CountryData has a built-in projection conversion for its country polygon:     *)
(* CountryData["USA", {"Polygon","projection"}]                                  *)
(* but it results in a grid system different from the GeoGridPosition conversion *)
(* So, I'll perform the transformation in the same way as above                  *)
USAplot = 
 Graphics[{FaceForm[LightGray], EdgeForm[Darker@Gray], 
   Apply[
      GeoGridPosition[GeoPosition[{#2, #1, 0}],"Mercator"][[1, ;; 2]] &,     
      CountryData["USA", "Polygon"], 
      {3}
   ]}]

cityPlot = BubbleChart[lc, BubbleSizes -> {0.01, 0.1}];
Show[USAplot, cityPlot]

You can prevent this from happening by borrowing the background graphics' options when making the BubbleChart plot:

cityPlot = 
 BubbleChart[lc, AbsoluteOptions[USAplot, {AspectRatio, PlotRange}],BubbleSizes -> {0.01, 0.1}];

Show[USAplot, cityPlot, AbsoluteOptions[USAplot, {AspectRatio, PlotRange}]]

Many legends can be easily created using combination of Row and Column, where for example, each Row contains two elements: a Graphics expression, and a text label which need not be wrapped in Text.

However, since the label you want involves size coding, this method doesn't work. Further, you want to match the size of Disks in the main graphics to the size in the legend, so I recommend using Epilog or Inset.

The basic layout of the legend itself is straightforward, up to scaling:

 Graphics[{Lighter@Blue, Table[{Disk[{i, 0}, Sqrt[i/Pi]/3],
    Black, Text[Style[i, FontSize -> 14], {i, -2/3}]}, {i, 1, 5}]}]

The scaling factor is up to you. For area-coding, the Sqrt is necessary since Disk is determined by radius, but typically the factor Pi is irrelevant since that's a constant parameter that only affects absolute size not relative size. Note the additional factor of 3 shrinking to avoid the disks From overlapping.

Here's a more compact legend using stacked Circles:

Graphics[{Lighter@Blue, Table[{Circle[{0, Sqrt[i/Pi]}, Sqrt[i/Pi]],
    Lighter@Gray,
    Line[{{1/4, 2 Sqrt[i/Pi]}, {2 - 1/4, 2 Sqrt[i/Pi]}}],
    Black,
    Text[Style[i, FontSize -> 14], {2, 2 Sqrt[i/Pi]}]}, {i, 1, 10}]},
 ImageSize -> 300
 ]