We consider an inverse scattering problems for the biharmonic Schrödinger operator Δ2 + A · ∇ + V in three dimensions. By the Helmholtz decomposition, we take A = ∇p + ∇ ×ψ. The main contributions of this work are twofold. First, we derive a stability estimate of determining the divergence-free part ∇ ×ψ of A by far-field data at multiple wavenumbers. As a consequence, we further derive a quantitative stability estimate of determining −12∇⋅A+V. Both the stability estimates improve as the upper bound of the wavenumber increases, which exhibit the phenomenon of increased stability. Second, we obtain the uniqueness of recovering both A and V by partial far-field data. The analysis employs scattering theory to obtain an analytic domain and an upper bound for the resolvent of the fourth order elliptic operator. Notice that due to an obstruction to uniqueness, the corresponding results do not hold in general for the Laplacian, i.e., Δ + A · ∇ + V. This can be explained by the fact that the resolvent of the biharmonic operator enjoys a faster decay estimate with respect to the wavenumber compared with the Laplacian.