How do I fix this heart?

The problem is simple. If we look at that horseshoe region, we get imaginary numbers. For part of it, it should belong to our heart. In that region, if we were to evaluate our function (by math, not by programming), the imaginary parts of the function cancel. So it should look like this (generated in Mathematica):

Basically, the function for that part is almost identical; we just have to do arithmetic with complex numbers instead of real numbers. Here's a function that does exactly that:

private static double topOfHeart(double x, double y, double temp, double temp1) {
    //complex arithmetic; each double[] is a single number
    double[] temp3 = cbrt_complex(temp, sqrt(-temp1));
    double[] part1 = polar_reciprocal(temp3);
    part1[0] *= 36 * cbrt(2) * pow(y, 3);
    double[] part2 = temp3;
    part2[0] /= (10935 * cbrt(2));
    toRect(part1, part2);
    double looseparts = 4.0 / 9 - 4.0 / 9 * pow(x, 2) - 4.0 / 9 * pow(y, 2);
    double real_part = looseparts + part1[0] + part2[0];
    double imag_part = part1[1] + part2[1];
    double[] result = sqrt_complex(real_part, imag_part);
    toRect(result);

    // theoretically, result[1] == 0 should work, but floating point says otherwise
    if (Math.abs(result[1]) < 1e-5) {
        return result[0];
    }
    return -1;
}

/**
 * returns a specific cuberoot of this complex number, in polar form
 */
public static double[] cbrt_complex(double a, double b) {
    double r = Math.hypot(a, b);
    double theta = Math.atan2(b, a);
    double cbrt_r = cbrt(r);
    double cbrt_theta = 1.0 / 3 * (2 * PI * Math.floor((PI - theta) / (2 * PI)) + theta);
    return new double[]{cbrt_r, cbrt_theta};
}

/**
 * returns a specific squareroot of this complex number, in polar form
 */
public static double[] sqrt_complex(double a, double b) {
    double r = Math.hypot(a, b);
    double theta = Math.atan2(b, a);
    double sqrt_r = Math.sqrt(r);
    double sqrt_theta = 1.0 / 2 * (2 * PI * Math.floor((PI - theta) / (2 * PI)) + theta);
    return new double[]{sqrt_r, sqrt_theta};
}

public static double[] polar_reciprocal(double[] polar) {
    return new double[]{1 / polar[0], -polar[1]};
}

public static void toRect(double[]... polars) {
    for (double[] polar: polars) {
        double a = Math.cos(polar[1]) * polar[0];
        double b = Math.sin(polar[1]) * polar[0];
        polar[0] = a;
        polar[1] = b;
    }
}

To join this with your program, simply change your function to reflect this:

if (temp1 < 0) {
    return topOfHeart(x, y, temp, temp1);
}

And running it, we get the desired result:

It should be pretty clear that this new function implements exactly the same formula. But how does each part work?

double[] temp3 = cbrt_complex(temp, sqrt(-temp1));

cbrt_complex takes a complex number in the form of a + b i. That's why the second argument is simply sqrt(-temp1) (notice that temp1 < 0, so I use - instead of Math.abs; Math.abs is probably a better idea). cbrt_complex returns the cube root of the complex number, in polar form: r eiθ. We can see from wolframalpha that with positive r and θ, we can write an n-th root of a complex numbers as follows:

And that's exactly how the code for the cbrt_complex and sqrt_complex work. Note that both take a complex number in rectangular coordinates (a + b i) and return a complex number in polar coordinates (r eiθ)

double[] part1 = polar_reciprocal(temp3);

It is easier to take the reciprocal of a polar complex number than a rectangular complex number. If we have r eiθ, its reciprocal (this follows standard power rules, luckily) is simply 1/r e-iθ. This is actually why we are staying in polar form; polar form makes multiplication-type operations easier, and addition type operations harder, while rectangular form does the opposite.

Notice that if we have a polar complex number r eiθ and we want to multiply by a real number d, the answer is as simple as d r eiθ.

The toRect function does exactly what it seems like it does: it converts polar coordinate complex numbers to rectangular coordinate complex numbers.

You may have noticed that the if statement doesn't check that there is no imaginary part, but only if the imaginary part is really small. This is because we are using floating point numbers, so checking result[1] == 0 will likely fail.

And there you are! Notice that we could actually implement the entire heart function with this complex number arithmetic, but it's probably faster to avoid this.