improving speed of this numpy-based diffraction calculator
I managed to make your script a bit faster. See the comments in the
source code:

import numpy as np
import matplotlib.pyplot as plt
import time

# Makes it so numpy always outputs the same random numbers which is
# useful during development. Comment out the line in production.
np.random.seed(1234)

def main():
    twopi, img = 2*np.pi, np.complex(0, 1)

    wavelength = 0.08
    k0  = twopi/wavelength

    z_noise = 0.05 * wavelength

    N, M = 100, 50
    x = np.arange(-N, N+1)

    # Use random Z noise for now
    Z = z_noise * np.random.normal(size = (x.shape[0], x.shape[0]))

    A = np.linspace(0, 0.2, M)

    tstart = time.time()
    answers = []

    # X and Y with shapes (x, 1) and (1, x) contains the same data.
    X = x[np.newaxis,:]
    Y = x[:,np.newaxis]

    # Compute squared Euclidean distances with shape (x, x).
    Rsq = X**2 + Y**2

    # Compute other distances with shape (M, M).
    a = A[:,np.newaxis]
    b = A[np.newaxis,:]
    SQ = 1 + np.sqrt(1 - a**2 - b**2)

    # Add noise thing (M, M, x, x).
    SQ_Z = SQ[:,:,np.newaxis,np.newaxis] * Z

    # Gets shape (M, 1, x) and (M, x, 1).
    A_X = A[:,np.newaxis,np.newaxis] * X
    A_Y = A[:,np.newaxis,np.newaxis] * Y

    # Calculates A*X + A*Y with shape (M, M, x, x).
    A_X_Y = A_X[:,np.newaxis] + A_Y[np.newaxis,:]

    # e^{-j*phases} with shape (M, M, x, x).
    A_stuff = np.exp(-img * k0 * (A_X_Y + SQ_Z))
    for w_beam in (1, 2, 4, 8):
        coeff = np.exp(-Rsq/w_beam**2)
        # Sums over the last two dimensions, shape becomes (M, M).
        E = np.sum(A_stuff * coeff, axis = (2, 3)) / w_beam**2
        answers.append(E)

    print(f'Time: {time.time() - tstart:.3}')
    # Print checksum. Useful for debugging.
    print(sum(np.sum(E) for e in answers))
    if True:
        plt.figure()
        for i, E in enumerate(answers):
            plt.subplot(2, 2, i+1)
            plt.imshow(np.log10(np.abs(E)), vmin=0.0001)
            plt.colorbar()
        plt.show()

main()


I have wrapped all your code in a main function. It is good practice
to do so, even for small scripts, but it makes it easy to see if
there's any dependencies on global variables.

The script runs about four times faster. Most of the speedup comes
from replacing your explicit loop with vector functions. So instead of

for i, a in enumerate(A):
    EE = []
    for b in A:
        e = E_ab(a, b, positions, w_beam, k0)
        EE.append(e)
    E.append(EE)


those lines now read

coeff = np.exp(-Rsq/w_beam**2)
E = np.sum(A_stuff * coeff, axis = (2, 3)) / w_beam**2
answers.append(E)


To get the data in a format that vector functions can handle I've
added a lot of
broadcasting. It
can be used to compute, for example, the cartesian
product
of arbitrary numpy arrays.

I think someone who understands the math (which I don't) could
optimize your script even more. Some of the calculations seem
redundant to me, but I have no idea how to factor them out.