I am attempting to find a mathematical curve $f(x,y)$ that describes approximately the contour plot of a bimodal constructed from two Gaussian distributions centered at $(0, \pm y_0)$ and $\sigma = 1$.
The curve $f(x,y)$ should resemble an ellipsis for points near and far away from $(0, \pm y_0)$ and transform into a type of Lemniscate for points in between.
As an example, see the following picture with $y_0 = 1.3$: 
Code to generate the image:
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
#Draw the bimodal distribution
N = int(1e8)
y_0 = 1.3
rng = np.random.default_rng(1)
x_g = rng.normal(size=2*N)
y_g = np.concatenate([rng.normal(loc = y_0, size=N), rng.normal(loc = -y_0, size=N)])
#Bins
x_max = 3
Bins_x = np.linspace(-x_max, x_max, 100)
Bins_y = np.linspace(-x_max, x_max, 100)
binCenter_x = (Bins_x[1:]+Bins_x[:-1])/2
binCenter_y = (Bins_y[1:]+Bins_y[:-1])/2
#Construct histogram
z_counts, edges = np.histogramdd((x_g, y_g), bins = [Bins_x, Bins_y])
#Plot
im = plt.pcolormesh(binCenter_x, binCenter_y, z_counts.T, cmap = 'plasma_r', zorder = -1,rasterized=True)#, norm = LogNorm(vmin=1, vmax=np.max(z_counts)))
plt.contour(binCenter_x, binCenter_y, gaussian_filter(z_counts, sigma=0.5).T, levels=10)
plt.colorbar(im, label = 'Counts')
plt.gca().set_aspect(1./plt.gca().get_data_ratio())
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('bimodal.png', dpi = 300)
The curves you are describing are conceptually identical to the Cassinian oval, as described in A Catalog of Special Plane Curves, J. Dennis Lawrence, Dover, 1972, pp. 153-154. A description can also be found here.
The Cartesian equation of these curves is given by
$$ (x^2+y^2+a^2)^2=b^4+4a^2x^2 $$
and the polar equation is
$$ r^4-2a^2r^2\cos2\theta=b^4-a^4 $$
These curves have the property that if $a=b$ the oval is lemniscate of Bernoulli. If $a<b$, there is one loop, and if $a\ge b$, two loops. The only difference between your curves and these is a $90^{\circ}$ rotation.