How do I apply BDF2 in a STRANG splitting

I have a 3D diffusion equation that I want to solve using a time splitting (2D+1D). Assume that $A$ is the 2D discrete diffusion operator and $B$ is the 1D discrete diffusion operator.

I want to use a STRANG splitting to ensure a global 2nd order of precision :
begin{equation}
frac{partial f}{partial t}=frac{1}{2}Af
frac{partial f}{partial t}=Bf
frac{partial f}{partial t}=frac{1}{2}Af
end{equation}
And for each time integration, I want to apply the BDF2 implicit formula to ensure stability and a local 2nd order precision:
begin{equation}
frac{partial f}{partial t}approx frac{f^{n+1}-frac{4}{3}f^n+frac{1}{3}f^{n-1}}{frac{2}{3}Delta t}
end{equation}
My question is : for the time iteration $nDelta t$, what is the $f^{n-1}$ of each stage of the splitting ? Is it the state of $f$ at the previous time iteration or the previous splitting stage ?

To be more literal, do I follow this algorithm ?
begin{equation}
f0,f1rightarrow(frac{Delta t}{2} A) f21
f1,f21 rightarrow(Delta t B) f22
f21,f22rightarrow(frac{Delta t}{2} A) f2
end{equation}
and for the following iteration we initialize $f0=f22$ and $f1=f2$.

Or this one ?
begin{equation}
f0,f1rightarrow(frac{Delta t}{2}A) f21
f0,f21 rightarrow(Delta t B) f22
f0,f22rightarrow(frac{Delta t}{2} A) f2
end{equation}
and for the following iteration we initialize $f0=f1$ and $f1=f2$.

Is there another way to ensure an unconditionnaly stable time integration splitting with a global and local 2nd order ?