Just wondering if anyone knows the history behind this discovery? The BCS mean field theory suggests a breaking of our $U(1)$ symmetry; hence, a complex order parameter makes sense in that regard. How was it first realized however?
I don't have a definitive historical answer, but here are some pointers.
The idea of using order parameters to describe symmetry-breaking phase transitions (in particular second-order or continuous phase transitions) is due to Landau and dates back to 1937. His central observation was that a symmetry is either broken or not, and that if the phase transition is continuous (no jump), then the (potential) energy near the phase transition must be able to be described as an expansion in the order parameter, which must follow the properties set out by the broken symmetry.
This was expanded by Ginzburg-Landau in 1950 to include deviations/fluctuations in the order parameter from a perfectly homogeneous average value. Not only did they include a gradient term for the order parameter, but also the vector potential by minimal coupling/Peierls substitution, to handle the magnetic field, as they were explicitly interested in describing superconductivity. The Ginzburg-Landau theory, even though it is a mean-field theory, is surprisingly accurate in many accounts for many superconductors. So probably the best answer to your question is the Ginzburg-Landau paper.
The concept of off-diagonal long-range order (ODLRO) was introduced by O. Penrose in 1951 and O. Penrose & Onsager in 1956. In modern terms, this is giving the two-point correlation function as $\langle O^\dagger(x) O(x') \rangle$, where $O(x)$ is the order parameter, and in this form applies to any long-range order due to symmetry breaking. ODLRO was important in establishing superfluidity as a quantum-coherent phenomenon (Bose-Einstein condensation), and as a symmetry- breaking phase transition.
Josephson predicted the effects named after him in 1962, because he was convinced the "$F$ function", or the order parameter, was a real, physical field, and not simply a mathematical tool. The Josephson effect is the best experimental evidence for taking the $U(1)$ order parameter seriously.
Nambu developed his theory of broken symmetry, culminating in what are now known as Nambu-Goldstone bosons, explicitly based on superconductivity, in 1960-1. This led to the understanding of particle physics in terms of broken symmetries, although the concept of order parameter is less important in that field. But much of the mathematical machinery, in particular relating to gauge fields, follows from this.
edit: I left out the BCS paper initially. I went through it again, and indeed it does not talk at all about broken symmetry. It does not reference Ginzburg-Landau directly, nor ODLRO. Of course, the BCS wavefunction is the "macroscopic wavefunction" in the Landau sense, but it is not discussed in any detail. The only thing related to broken symmetry is the phrase:
If for the moment we relax the requirement that the wave function describes a system with a fixed number of particles
Furthermore, Bardeen did at first not believe the Josephson effect could be real, so this way of looking at it was definitely not on his mind. The term order parameter is mentioned once, in the Conclusions, but this is only about the amplitude (the gap function), which says nothing about the $U(1)$ symmetry breaking.