ψ[x_, y_] = 
      Exp[I (x kx + y ky) - 
          ((x - rx)^2 + (y - ry)^2)/(4 ϵ^2)]

P[x_, y_] = Conjugate[ψ[x, y]] ψ[x, y] // ComplexExpand

Integrate[P[x, y]  /.  rx -> 0  /.  ry -> 0, x, y]

Limit[Erf[x] , x -> Infinity] - Limit[Erf[x], x -> -Infinity]

Integrate[
      P[x, y], {x,  -Infinity,  Infinity}, {y,  -Infinity,  
        Infinity}]

<< Calculus`FourierTransform`

ψx[x_] = ψ[x, y]  /.  y -> ry /. ry -> 0

ψy[y_] = 
      ψ[x, y]  /.  x -> rx /. rx -> 0 // ComplexExpand

Integrate[
      ψx[x] Exp[2  Pi I f x], {x, -Infinity, Infinity}, ]

InverseFourierTransform[Out[50] // ComplexExpand, q, x]

<<Statistics`MultinormalDistribution`

\(PDF[
          MultinormalDistribution[{rx, 
              ry}, {{ϵ^2, 0}, {0, ϵ^2}}], {x, y}]^
        \((1/2)\)
 // PowerExpand\)

StandardDeviation[
      MultinormalDistribution[{rx, 
          ry}, {{ϵ^2, 0}, {0, ϵ^2}}]]