Are there experiments that could show that light waves resemble more say square waves than sine waves? Are there experiments that could show that light waves resemble more say square waves than sine waves?
Temporarily looking at your example of a square wave, a square wave of spatial wavenumber $k$ can be represented in a Fourier expansion as a superposition of sine waves with spatial wavenumbers $k,3k,5k...$ etc. Assuming that we are doing our measurements in a medium which has linear response to the fields (which is satisfied in most materials for fields with power densities $\ll10^8\mathrm{V/m}$, which virtually all optics experiments satisfy), the superposition principle holds, and so without loss of generality we can examine the effect of an experimental apparatus on each individual basis function, and then add together the results.
It's known that individual sine waves will diffract from a diffraction grating at well-defined angles, which is how monochromators work. So, if we had a square-wave source, we would see multiple harmonics split by the grating. Alternately, we could use a prism. In either case, if (quasi)monochromatic light really was a square wave, you'd be able to see the harmonics, but you don't.
(From a comment): You say "It's known that individual sine waves will diffract from a diffraction grating at well-defined angles", but how do we know that? Individual waves will diffract from a diffraction grating at well-defined angles, but how do we know these are sine waves?
It's not quite correct to say that "Individual waves will diffract from a diffraction grating at well-defined angles" in the sense that the diffraction pattern will take the form of a "picket fence" of single peaks. Read any intro optics text and you'll see examples (such as multiple-slit diffraction and diffraction gratings) where they temporarily assume that the electric field takes the form of sinusoids, and then go on to show that sinusoid waves diffract at integer multiples of an angle. The reasoning does not necessarily apply to arbitrary waves (but due to linearity, arbitrary waves can be handled by Fourier decomposition).
(From a comment): Also couldn't we represent a sine wave as an infinite superposition of square waves as well? Or at least as a superposition of periodic functions which are not sine nor cosine? In which case by your reasoning we would be able to see the harmonics even with a sine wave
While it's tempting to use verbal reasoning to come to conclusions in physics, you really have to do the math to make sure your words are correct. You are correct in saying that a decomposition in terms of sines is not necessary; by linearity, we can actually compute the physics using any periodic basis we want, such as square waves, and because of the linearity, the end result will always be the same.
To illustrate this, let's see what exactly happens when we attempt to diffract a sine wave and a square wave on a multiple-slit setup with 61 slits. Let the $k$th slit be located at a height $kd$ where $d$ is the slit spacing, and let the screen be located at a distance $R$ from the slits. At a location of height $y$ on the screen, the distance between the $k$th slit and the location is
$$d(k,y)=\sqrt{R^2+(y-kd)^2}\approx \sqrt{d^2 k^2+R^2}-\frac{d k y}{\sqrt{d^2 k^2+R^2}}$$
and for an incident sine wave the field amplitude hitting the point becomes
$$E(k,y)\propto\exp\left(\frac{2\pi i d(k,y)}{\lambda}\right)$$
and so the light intensity at that point becomes
$$I(y)=\left|\sum_{k=-30}^{30}E(k,y)\right|^2$$
which basically looks like this:
expr = Sum[
Exp[2 \[Pi] I (Sqrt[d^2 k^2 + R^2] - (d k y)/Sqrt[
d^2 k^2 + R^2])/\[Lambda]], {k, -30, 30}];
Plot[Abs[expr /. {R -> 1, \[Lambda] -> 0.0000025,
d -> 0.00002}], {y, -0.3, 0.3}, PlotRange -> All,
PlotPoints -> 60, ImageSize -> 900, AspectRatio -> 0.3]
Meanwhile, we can do the exact same thing with a square wave beam hitting the slits. We have
$$E(k,y)\propto\mathrm{SquareWave}\left(\frac{d(k,y)}{\lambda}\right)$$
and again the light intensity at that point becomes
$$I(y)=\left|\sum_{k=-30}^{30}E(k,y)\right|^2$$
which basically looks like this:
expr = Sum[
SquareWave[(Sqrt[d^2 k^2 + R^2] - (d k y)/Sqrt[
d^2 k^2 + R^2])/\[Lambda]], {k, -30, 30}];
ListLinePlot[
Table[Abs[
expr /. {R -> 1, \[Lambda] -> 0.0000025, d -> 0.00002}], {y, -0.3,
0.3, 0.00005}], PlotRange -> All, AspectRatio -> 0.3,
ImageSize -> 900]
Notice that in addition to the main peaks, there are smaller sideband peaks with spacings which have spacings which are factors of $1,3,5,7,...$ times smaller that of the fundamental sequence. These are the harmonics I mentioned earlier, which as I mentioned before can be derived by considering the decomposition of the square wave into a sinusoid basis.
However, note that nowhere in the code above used to generate the image did I use sine waves! It was entirely based on a direct square wave summation. This is a good illustration of the fact that the choice of basis you choose to represent linear physics in is irrelevant.
