Advice on designing a digital filter that doesn't have phase-sensitive edge artifacts?
Your basic problem is that filtfilt (and most other linear filtering routines) take filters that are designed for infinitely long time expanses, and apply them to a chunk of data as if the data were extended infinitely in both directions with zeros.
So you have a legitimate bandpass filter, and it's &quot;seeing&quot; a legitimate jump in the signal at the starting point of your signal.
There's three basic approaches you can take; the first two are ad-hoc and easy, the third one is difficult if you're starting from first principles.  It's certainly been solved out there someplace, but a brief search here on &quot;filter finite-length data&quot; didn't find me joy.
Approach 1: window the input data
Take your input data, and multiply it by something that'll make it taper off at the ends.  I.e. a ramp from 0 to 1 over 10 samples at each end, or $\frac{1}{2}\left (1 - \cos \frac{\pi n}{N} \right)$ for N samples at each end (suitably reversed on the left end).  You'll have some artifacts (a rising sine wave isn't the same as a steady one, after all), but they'll be attenuated.
Here's python code implementing a cosine edge attenuation with the ability to customize what central % of the signal is kept as 1:
def attenuate_edges(signal,time,edge_attenuation_percent):
  start = int(np.floor(len(time)*edge_attenuation_percent))
  end = int(len(time)-start)
  ramp = (1-np.cos(np.pi*(np.arange(start)/start)))/2
  edge_attenuator = np.ones(len(time))
  edge_attenuator[0:start] = ramp
  edge_attenuator[end:len(time)] = np.flip(ramp)
  return(signal*edge_attenuator)

Approach 2: Trim the output data
Do what you're doing now, and lop off the nastiness at the ends.  This is probably the easiest, and if you can just collect a bit more data, doesn't lose you anything.
Approach 3: Do a proper estimate of the interfering signal, and subtract it out
This will be fun-fun if you love math and have the time.  Basically you'll use the fact that the value of your interfering signal at time $n$ correlates in a specific way with the values of your interfering signal at time $k$ for all values of $n$ and $k$ in your data set.  You'll probably end up with something that looks a lot like a Wiener or a Kalman filter, that takes the end-effects into account.  Your estimate will be worse at the ends, but this will show up as a bit of noise on the ends -- not as those honkin' big pulses.
If I couldn't figure out the search terms for this, it would take me a day to do and another day to verify, and supposedly an expert.  OTOH, Gauss or Laplace probably invented it in the 19th century, and may even have thought it important enough to write down, somewhere.  So I'm sure the method exists.