Is jaggedness truly an "aliasing" effect?

The word aliasing has a precise mathematical meaning. Let a function $f:[-1,1]\to\mathbb{R}$ be given. Now suppose you want to approximate $f$ as best as possible wrt the standard complex inner product on $[-1,1]$: $\langle g, f \rangle = \int_{-1}^{1} \overline{g(x)} f(x)\,dx$. If you have an orthogonal Schauder basis $(b_k)_{k\in\mathbb{Z}}$ for $L_2([-1,1])$, where $b_k:[-1,1]\to \mathbb{C}$, you can approximate any "function" in $L_2$ as a linear combination of the basis functions: $\exists (c_k)_{k\in\mathbb{Z}} \,:\,\|f - \sum_{k\in\mathbb{Z}} c_k b_k\|_2 = 0$. Since the basis is orthogonal it holds that $$c_k(f) = \frac{\langle b_k, f\rangle}{\langle b_k, b_k\rangle} = \frac{\int_{-1}^{1} \overline{b_k(x)}f(x)\,dx}{\int_{-1}^{1} |b_k(x)|^2\,dx}.$$

One such representative is the trigonometric basis $b_k(x)=\exp(i \pi k x)$. One can easily verify that it is orthogonal wrt the above inner product and that $\|b_k\|^2 = 2$. Then the series $$S(f)(x) = \sum_{k\in\mathbb{Z}} c_k(f) b_k(x)$$ is known as the Fourier series of $f$. The first issue with the above is that we may be unable to compute the integrals for the coefficients even if we truncate the sum. You could then choose to approximate the integrand by a piecewise constant function of $2N$ pieces - then the integral of the latter is just the sum of the areas of the rectangles: $$c_k = \frac{1}{2} \int_{-1}^1\exp(-i\pi k x)f(x)\,dx \approx \frac{1}{2}\sum_{j=-N}^{N-1} \exp(-i\pi k j/N) f(j/N) \frac{2}{2N} = \tilde{c}_k.$$ We can try to express the approximated coefficients $\tilde{c}_k$ by the true coefficients $c_{\ell}$. To that end suppose $S(f) = f$ and use the series representation of $f(j/N)$: \begin{align} \tilde{c}_k &= \frac{1}{2N}\sum_{j=-N}^{N-1} \exp(-i\pi k j/N) f(j/N) \\ &= \frac{1}{2N} \sum_{j=-N}^{N-1} \exp(-i\pi kj/N) \sum_{\ell=-m}^m c_{\ell} \exp(i\pi \ell j/N) \\ &= \frac{1}{2N}\sum_{\ell=-m}^m c_{\ell} \sum_{j=-N}^{N-1} \exp(i\pi (\ell-k) j/N) \\ &= \frac{2N}{2N} \left(\sum_{r\in\mathbb{Z}} c_{2rN+k}\right). \end{align} The last equality holds because if $\ell-k=2rN$ (is a multiple of $2N$) then $$\exp(i \pi (k+\ell) j / N)=\exp(i \pi 2rN j/N)=\exp(i2\pi rj) = 1.$$ And we have 2N summands, so the sum is $2N$. On the other hand when $p=k+\ell\ne 2rN$ the sums are zero, this is the same discrete orthogonality principle used in the DFT \begin{align} \sum_{j=-N}^{N-1} \exp(i \pi p j/N) &= \sum_{j=0}^{2N-1} \exp(i \pi p (j-N)/N) \\ &= \exp(-i\pi)^p\sum_{j=0}^{2N-1} \exp(i \pi p/N)^j \\ \\ &= (-1)^p \frac{\exp(i\pi p/N)^{2N} - 1}{\exp(i\pi p/N)-1} = 0. \end{align} The last holds from the formula for the sum of geometric series.

Ultimately the key observation is that $\tilde{c}_{k} = \sum_{r\in\mathbb{Z}} c_{2rN+k}$. That is, with this sampling we cannot distinguish between waves of frequency $k$ and $2rN+k$. This is precisely what aliasing refers to. The first frequency at which this occurs is $k=N$ since we cannot distinguish it from $k-2N = N-2N = -N$. This is precisely the Nyquist limit as we had $2N$ samples and issues occur at $N$ which is half of the sampling rate $2N$. This also means that if our signal is bandlimited: $c_k=0$ for $|k|\geq N$ we have no issues ($k\geq N$ and $k<-N$ is also fine in the discrete setting with the above formulation).

TLDR

The takeaway of this lengthy example is that aliasing is the result of numerical integration error in your best approximation formulation. For example if you have to draw a triangle, ideally you would color a pixel based on the area of the triangle covering it, but that is expensive, so modern GPUs instead just sample a few locations to approximate coverage. This causes jaggies. Similarly if we have a high resolution texture you can see aliasing clearly. That's why you have mipmaps. Analytic filtering/integration of procedural textures is sometimes also possible.

Finally, if you have a very low resolution screen then even if you were to integrate perfectly within each pixel you would still have a blocky image. The mathematical counterpart of this is truncating the series even if you were to exactly compute the coefficients. In practice, however, a lower screen resolution typically also implies a higher integration error as your rendering resolution is usually tied to the output resolution (unless you decide to render to a target with a resolution several times higher and then downsample, but then again edges have infinite frequencies). I also don't think screen resolution is an issue nowadays since at 4K or 8K I woulad argue your eyes' PSF already filters out any frequencies that would cause aliasing, so the only thing causing aliasing there is integration error from rendering with an insufficient sample rate or/and insufficient prefiltering.

So in theory the answer to your question is it depends. But in practice jaggedness is almost always actual aliasing due to integration error. Regarding different reconstruction functions - yes they will produce different artifacts, the Fourier basis has actually quite unpleasant and global oscillatory artifacts (see Gibbs' phenomenon), on the other hand highly localized and regular stuff like box functions tends to result in more jaggies as you noted. Ideally one would pick something inbetween like Lanczos or w/e, but higher order operators typically suffer from oscillations and halos just like the Fourier basis. Ultimately it's a tradeoff between what kinds of artifacts you want your integration error (or interpolation assumptions) to show up as - blur vs sharpness against halos/oscillations. The sinc function is actually a terrible choice in practice - windowed/weighted sinc functions are more acceptable. Also note that the standard shannon theorem applies only to regular sampling. One can do much better with adaptive sampling and appropriate functions such as wavelets.