When solving 1d particle in a box, the probability density is said to be proportional to $|\psi|$, but when solving 3d orbitals, the probability density is said to be proportional to $|\psi|^2 r^2$. Why this difference?
[Physics] Why wavefunction is sometimes multiplied by the radius to get probability density
orbitalsprobabilitywavefunction
Related Solutions
All the definitions you've posted are correct, and they aren't in conflict with each other, although they are a bit imprecise. I'll try to explain in more detail what these concepts are.
Probability
I'll assume we don't need to attempt to define what probability is. :) I'll just note that it's formalized in mathematics under the aegis of measure theory.
Probability Density
Probability density is a concept that naturally arises whenever you talk about probability in connection with a continuous variable, such as position of a particle. In contrast to discrete probability (such as idealized coin-flipping or dice-rolling), we can't directly assign a probability to each individual outcome. Probabilities have to sum to one, but there are infinitely many possible positions for a particle, so any individual position would have to have probability zero. So how can we talk about probability spread over a continuous domain?
Well, we can talk about the probability of finding the particle within some region. If the region covers a nonzero length/area/volume/whatever, we can expect it to have a nonzero probability. But the exact probability will generally depend on the size/shape/position of the region in a complicated way.
Suppose we consider a very small region. To be concrete let's suppose we're talking about a 1D position and our small region is some interval $[x, x+dx]$. The interval has width $dx$. It stands to reason that a small interval will have a small probability. Moreover, if probability is "spread smoothly" in a sense that I won't bother making precise (hey, this is Physics.SE, not Math.SE!), then for small intervals the probability will be proportional to the width of the interval. If you halve the size of the interval you get half as much probability. Therefore $$dP = p(x) \, dx$$ where $dP$ is the small probability of finding the particle in this small region, and $p(x)$ is some proportionality constant. I've written it as a function of $x$ because the proportionality could vary from point to point.
Well, $p(x)$ is called the probability density, and it has the units of probability per unit length. (Or unit area, or volume, depending on how many dimensions we do this in.) It's a lot like mass density, which is mass per unit length/area/volume, or charge density is charge per unit length/area/volume, etc...
So, now you can imagine calculating the probability to find the particle in any arbitrary region, by slicing up that region into many small ones and summing up $dP$, which is $p(x) \, dx$, for all of them. Well, this is just integrating: $$P(\text{particle is in region $R$}) = \int_R p(x) \, dx$$ In particular, the particle has to be somewhere, so $p(x)$ must integrate to one: $$P(\text{particle is somewhere}) = 1 = \int_{-\infty}^\infty p(x) \, dx$$ But $p(x)$ itself isn't limited to 1. It isn't a probability but a probability per unit length/area/volume, so $p(x) > 1$ just means a region where probability is concentrated. (Much like you can drive your car at 100 km/hr without actually driving 100 km.)
The Wavefunction
All of the above was just standard probability theory, nothing to do with quantum mechanics. Now, in QM we try to predict the probability density for a particle's position (or momentum, or energy, or whatever). We could try to do this by writing an equation for how $p(x)$ changes over time, but it turns out that doesn't give us enough information; there are situations where particles start with identical $p(x)$ but do different things as time goes on.
It's found that we do get enough information to make predictions if we write an equation for a complex-valued function $\psi(x)$, and derive the probability density from it as $$p(x) = \psi^*(x) \psi(x)$$ The way the complex phase of $\psi(x)$ varies from point to point encodes additional information about the particle's momentum, which is necessary to predict its future behavior. It has units of the square root of a probability density, which is a bit weird but perfectly mathematically acceptable. This is of course the wavefunction, and the equation that determines how it varies is the Schrödinger equation.
No, you can't.
The function $\psi\in\mathbb C$ has two real degrees of freedom; they are coupled and dynamical (non-gauge). On the other hand, the function $\rho\in\mathbb R$ has one real degree of freedom. It is impossible to reduce the dynamics of the system from two variables to one variable without losing information in the process.
(But, in a formal sense: Yes, you can)
Let $\psi=\sqrt{\rho}\mathrm e^{iS}$, with $\rho,S$ a pair of real variables. You may write the Schrödinger equation directly in terms of $\rho,S$ as (cf. Madelung or Bohm) \begin{equation} \begin{aligned} \frac{\partial\sqrt{\rho}}{\partial t}&=-\frac{1}{2m}\left(\sqrt{\rho}\nabla^2S+2\nabla\sqrt{\rho}\cdot\nabla S\right)\\ \frac{\partial S}{\partial t}&=-\left(\frac{|\nabla S|^2}{2m}+V-\frac{\hbar^2}{2m}\frac{\nabla^2\sqrt{\rho}}{\sqrt{\rho}}\right) \end{aligned} \end{equation}
As you can see, you cannot write an equation for $\rho$ alone, because its equation is coupled to a second unknown, $S$. Two real degrees of freedom, not one. Formally speaking, you may solve the equation for $S$ as a functional of $\rho$, and plug the result into the equation for $\rho$, thus obtaining an equation for $\rho$ alone. This is impractical because it is not really possible to solve for $S=S[\rho]$ in general terms, and even if we could, the functional would be highly non-local so the resulting equation for $\rho$ would be impossible to work with. The Schrödinger equation, written in terms of $\psi$, even if complicated, is as simple as it gets. Any other reformulation is way more cumbersome to use.
Best Answer
It's not "multiplied by $r^2$ to get the probility density". The issue is that the volume element in spherical coordinates is $$ \mathrm{d}V = r^2\sin(\theta)\mathrm{d}r\mathrm{d}\theta\mathrm{d}\phi$$ and since the probability to find a particle in a subspace $X\subset \mathbb{R}^3$ is $$ P(X) = \int_X\lvert \psi(r)\rvert^2\mathrm{d}V$$ by definition of a probability density, the quantity $r^2\lvert \psi(r)\rvert^2$ is what behaves like the "normal" probability density in flat coordinates: The probability to find the particle between $r_1$ and $r_2$ is proportional to $ \int_{r_1}^{r_2} r^2\lvert \psi \rvert^2\mathrm{d}r$.