We compute the electromagnetic form factors of the proton and neutron using lattice QCD with $N_f = 2 + 1 + 1$ twisted mass clover-improved fermions and quark masses tuned to their physical values. Three ensembles with lattice spacings of $a$=0.080 fm, 0.068 fm, and 0.057 fm, and approximately the same physical volume allow us to obtain the continuum limit directly at the physical pion mass. Several values of the source-sink time separation ranging from 0.5 fm to 1.5 fm are used, enabling a thorough analysis of excited state effects via multi-state fits. The disconnected contributions are analyzed using high statistics for the two-point functions combined with low-mode deflation and hierarchical probing for the fermion loop estimation. We study the momentum dependence of the form factors using the z-expansion and dipole Ansaetze, thereby enabling the extraction of the electric and magnetic radii, as well as the magnetic moments in the continuum limit, for which we provide preliminary results.