what is the 2D gravity potential?

You're right that if you take Newton's law of gravity as is and apply it to a 2D universe, you'll get an infinite result. So you do need to use a modified theory in two dimensions, or indeed in any number of dimensions other than three.

The proper way to do this is using general relativity, and if you apply GR to 2+1D spacetime, you get something that looks basically nothing like gravity as we know it. In particular, space is only distorted or "curved" where there is actually mass, unlike our universe where the distortion extends beyond the region that actually contains the mass. Because that distortion is what we recognize as gravity, in a 2+1D world there would be no gravitational attraction. The presence of mass would cause some geometrical oddities, but there would be no force acting between separated masses. For details, see e.g. "General relativity in a (2 + 1)-dimensional space-time".

Before GR was invented, on the other hand, physicists would have tried to generalize Newtonian gravity to other numbers of dimensions using Gauss's law for gravitation, which is exactly equivalent to Gauss's law for electric fields. As you may know, the electric version of Gauss's law states that the surface integral of electric field over a closed surface is proportional to the enclosed charge, so analogously, the surface integral of the gravitational field over a closed surface is proportional to the enclosed mass.

$$\int \vec{g}\cdot\mathrm{d}\vec{A} = -4\pi GM$$

In 3D, for a sphere enclosing a point mass, this gives you $(4\pi r^2\hat{r})\cdot(g\hat{r}) = -4\pi GM$, or $\vec{F}_g = m\vec{g} = -\frac{GMm}{r^2}\hat{r}$, which is recognizable as Newton's law of gravitation. In 2D space, you instead get $(2\pi r\hat{r})\cdot(g\hat{r}) = -4\pi GM$, which leads to

$$\vec{F}_{g\text{(2D)}} = m\vec{g} = -\frac{2GMm}{r}\hat{r}$$

So to do your integral with 2D Newtonian gravity, you would replace the exponent $\frac{3}{2}$ in the denominator of the integral with $1$, as well as changing the constant (but I'll skip that, since you've set it to 1 anyway).

$$F = \int_{-1}^1 \int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}} \tfrac{x+1}{(x+1)^2+y^2} \mathrm{d}y\, \mathrm{d}x = 2\pi$$

I may have done that integral wrong, but the point is, the answer is finite. You could also use Gauss's law directly to see it: a unit circle filled with material of unit surface density encloses a mass $\pi$, so the 2D Newtonian gravitational acceleration at the surface would be $-\frac{2GM}{r} = -2\pi G$.

From this modified force, you could determine that the potential energy would indeed have to have a logarithmic form, just using $U = -\int \vec{F}\cdot\mathrm{d}\vec{r}$.

The usual definition of gravity in different dimension makes it so that if space is 2d, then Newtonian gravity goes as 1/r not 1/r^2 (as David Zaslavsky points out, GR in 2+1 dimensions is even weirder--- gravity doesn't exert long range forces in 2d, but it creates point conical deficits. I will ignore GR in this answer), so that the potential is logarithmic. This is because the proper 2d 1/r gravity is harmonic on the punctured plane, just like 3d gravity. The orbits are not ellipses anymore, they aren't closed curves, just radially oscillations are with a different period than the angular oscillations. Also, all orbits are bound. This is not what you are asking about, but you should know it.

What you are doing is using the 3d potential for a mass density which has been squashed into a plane. The divergence you see is due to the fact that you have an infinite density on the plane--- you have a massive plate of thickness $\epsilon$ with a mass density per unit area, so an infinite mass density per unit volume.

Away from the circle edge of the plate, when you are close to the plate, it looks like an infinite planar mass density. In this case, the near-plate solution for the potential is that it looks like the absolute value function $|a|$ where a is the perpendicular distance from the infinite plate. This is the same as in electrostatics--- the solution for the gravitational field for an infinite plate with constant mass density is by Gauss's law. You have a gravitational field pointing toward the plate which is finite at the plate, but which has a finite jump at the plate, and the jump magnitude is proportional to the density per unit area. This means that the potential has a discontinuous first derivative, but it isn't infinite at the plate.

But right near the circular edge, the near plate solution is different. In the limit that you are very close to the edge, you are looking at a half-infinite plate with a finite mass density. This scaling limit can be solved using 2-d methods. If you put the plate in the z-x plane, and stretching along the negative x-axis, the potential is a function only of y and x, not of z, so there is no electric field along the z-direction. The electric field is a harmonic function away from the negative x axis, and it is a 2-d harmonic function. This means that G_x + iG_y is the real part of a holomorphic function (G is the gravitational field).

For a point mass, G_x + i G_y turns out to be the holomorphic (except for one point) function $1/z$. So each pole is the location of a mass. The finite mass density along the negative x-axis means that the harmonic function has a constant residue density along the negative x-axis. The residue density is the cut-discontinuity (as explained in this answer: Correlation function which has branch cut in momentum space ), so you find that the solution is the log function:

$$ E_x + i E_y = \log(z) $$

This is indeed divergent at z=0, but not on the plate. This explains why your integral is divergent, and tells you the asymptotic form of the divergence.

This divergence at the edge of a flat plate is not a surprise. It is even worse when you have a thin line, then the divergence is as 1/r. Plates have diverging fields, and thin needles have even more diverging fields. This is the principle of the lightning rod, you make a thin line at zero potential to cause strong electric fields which ionize and discharge the air. This is not quite the same, because your plate is not at equal potential, but it's the same general principle.

What is going on here? How is the gravity potential derived?

In the context of non-relativistic gravity, simply look to the fundamental solution for Laplace's equation. From the Wikipedia article "Fundamental solution":

Laplace equation

For the Laplace equation,

$[-\nabla^2] \Phi(\mathbf{x},\mathbf{x}') = \delta(\mathbf{x}-\mathbf{x}')$

the fundamental solutions in two and three dimensions are

$\Phi_{2D}(\mathbf{x},\mathbf{x}')= -\frac{1}{2\pi}\ln|\mathbf{x}-\mathbf{x}'|,\quad \Phi_{3D}(\mathbf{x},\mathbf{x}')= \frac{1}{4\pi|\mathbf{x}-\mathbf{x}'|}$

So yes, in 2D, the potential is logarithmic.

If so, why, and isn't it a paradox if point masses in 2D orbit according to a different law than coplanar point masses in 3D?

It's not a paradox; the force law between the two point masses is different in 2D from what it is in 3D. Consider a purely radial orbit. Would you consider it a paradox that colinear point masses in 3D don't follow a 1D force law?