I've just learned about moments of inertia in my physics I class, $I=\int{r^2dm}$. The math involving moments of inertia in relation to torque and angular momentum is clear to me. I'm looking for some intuition towards where the definition/idea of moments of inertia comes from? Where do we obtain that the moment of inertia of a point mass $I=mr^2$?
I'm interested in the logical development of these ideas.
knzhou has explicitly showed in his post how the form $\displaystyle\int r^2~\mathrm dm$ came to picture. So, I'll be not writing the same thing again. Nevertheless, since, OP is interested in the "logical development" of the idea, I deem it worthy to add a brief discussion on moment of inertia as a second-rank tensor in a more general case so that OP might get a view of the whole topic in short.
The angular momentum $\mathbf L$ of a rigid body with respect to a certain stationary point is expressed as:
$$\begin{align}\mathbf L &= \sum_i \mathbf r_i\times (m_i\mathbf v_i)\\ & =\sum_i \mathbf r_i \times m_i (\boldsymbol \omega \times \mathbf r_i)\end{align}$$ where $\mathbf r_i$ is the distance of the $i$ th particle from the fixed point.
Using the identity $\mathbf A\times \mathbf B\times \mathbf C= \mathbf B(\mathbf A\cdot \mathbf C)- \mathbf C(\mathbf A\cdot \mathbf B),$ we get $$\mathbf L = \sum_i \left[m_i~\boldsymbol \omega r_i^2 - m_i \mathbf r_i(\mathbf r_i\cdot \boldsymbol \omega)\right]\tag I$$
From $\rm (I)$ we get the $x$ component of $\mathbf L$ as
$$L_x = \sum_i\left[ m_i(r_i^2 -x_i^2)~\omega_x - m_ix_iy_i~\omega_y- m_ix_iz_i~\omega_z\right]$$
So, $L_z$ is a linear function each component of angular velocity.
Thus, $\mathbf L$ is linearly transformed from $\boldsymbol \omega\,.$
This can be expressed as $$\mathbf L = \mathsf I\boldsymbol\omega\tag{II}$$
where $\sf I$ is an operator which linearly transforms $\boldsymbol\omega$ to $\mathbf L\,.$ That is,
$$L_i = \sum_j ~I_{ij}~\omega_j$$
$$\therefore ~~~~~~~~~ \mathsf I \equiv \begin{bmatrix}I_{xx} & I_{xy}& I_{xz}\\ I_{yx}& I_{yy}& I_{yz}\\ I_{zx}& I_{zy}& I_{zz} \end{bmatrix} $$ where $$I_{ij} = \sum_ i m_i (r_i^2 ~\delta_{ij} - r_ir_j)\tag{III}$$ and $\mathbf L$ and $\boldsymbol\omega$ are column matrices.
From $\rm(III),$ it can be concluded that $I_{xy} = I_{yx};$ this means $\sf I$ is a symmetric matrix.
Any symmetric linear operator can be expressed in a unique diagonal form (the diagonal elements are unique but the order of the elements may be different) by orthogonal transformation of coordinates. The corresponding principle axes will be unique and orthogonal.
Now, kinetic energy $\textrm{KE}$ of the rigid body is $$\begin{align}\textrm{KE}&=\frac12 \sum_i m_i (\boldsymbol \omega \times\mathbf r_i)\cdot (\boldsymbol \omega \times\mathbf r_i)\\ &=\frac12\sum_i m_i \boldsymbol\omega\cdot (\mathbf r_i\times (\boldsymbol \omega \times\mathbf r_i))\\&= \frac12\boldsymbol\omega\cdot \mathbf L\\ &= \frac12 \boldsymbol\omega\cdot \mathsf I\boldsymbol \omega~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(\textrm{using}~\mathrm{II})\end{align}$$ which on expanding gives $$\begin{align}T &=\frac12( I_{xx}\omega_x^2 + I_{yy}\omega_y^2 + I_{zz}\omega_z^2) + \frac12\sum_i\sum_j(I_{ij}\omega_i\omega_j) \\ &= \frac12\sum_{i}I_{ii}~\omega_i^2 + (I_{xy}\omega_x\omega_y + I_{yz}\omega_y\omega_z + I_{zx}\omega_z\omega_x)~~~~~~~~~~~~~(I_{xy}= I_{yx}) \end{align}\,.$$
Now, when $x,y,z$ are the principle axes, then $$T=\frac12( I_{xx}\omega_x^2 + I_{yy}\omega_y^2 + I_{zz}\omega_z^2)\,. $$ The quadratic equation implies that the locus of angular velocities for which $\mathrm{KE} $ is constant must be an ellipsoid since energy is always positive.
This means all the off-diagonal elements of $\mathsf I$ will be zero and $\mathbf L \parallel\boldsymbol{\omega}\,.$
So, in principle axes, only the diagonal coefficients of moment of inertia are non-zero viz. $$L_i = I_{ii}\omega_i^2\,.$$
$\bullet$ Classical Mechanics by Goldstein.
Let's consider linear acceleration, governed by $F = ma$. If you have a system with many particles, and want to accelerate a system as a whole, the inertia depends on the total mass $$m_{\text{tot}} = \sum m_i = \int dm$$ where I just converted the sum over masses into an integral.
Now consider a single particle with mass $m$ rotating at radius $r$. Using $\tau = F r$ and $a = r\alpha$, the relationship between torque and angular acceleration is $$F = m a \quad \to \quad \frac{\tau}{r} = mr \alpha \quad \to \quad \tau = (mr^2) \alpha.$$ We define the moment of inertia to be $I = \tau/\alpha$, analogous to $m = F/a$. Then we see the moment of inertia of a particle is $mr^2$.
Now suppose you want to make a system rotate as a whole. Since each individual particle contributes a moment of inertia of $m_i r_i^2$, the total moment of inertia is just the sum, $$I_{\text{tot}} = \sum m_i r_i^2 = \int r^2 dm.$$ That's where the formula comes from.
