Simplest possible and most fundamental model of a solar cell comes from detailed balance of absorbed and emitted light.
This is not really a device model, which I think you are more interested in. But the electrical current that you get out of an ideal solar cell is the difference between the photon current the solar cell absorbes the photon current it emits,
$$ J(V) = q\left( N_a - N_e(V) \right) $$
Normally these are areal quantities: photons per square meter and current per square meter.
If you assume the appropriate blackbody spectrums, this is equivalent to the <PERSON>–Queisser model in the radiative limit (no non-radiative recombination).
Idealised model
We are going to calculate the current-voltage characteristic of an ideal Silicon solar cell.
Assuming the sun is a blackbody at temperature $T_s=5760K$, the photons per metre squares per second per Joule arriving at the Earth's surface is given by the Planck equation in terms of photon energy $\hbar\omega$,
$$ n_s(\hbar\omega) = f_s\left(\hbar\omega\right)^2\frac{2\pi}{c^2h^3}\frac{1}{\exp{\frac{\hbar\omega}{kT_s}}-1} $$
where the angular factor $f_s=6.8\times 10^{-5}$ considers the $1/r^2$ losses due to the Earth-Sun distance.
We are going to assume that the solar cell absorbs all sunlight above its bandgap energy $E_g$ and is completely transparent otherwise.
$$ a(\hbar\omega) = \begin{cases} 1,& \text{if } \hbar\omega \ge E_g\ 0, & \text{otherwise} \end{cases} $$
The number of photons absorbed from the sun is,
$$ N_a = \int_{0}^{\infty} a(\hbar\omega) n_s(\hbar\omega) d(\hbar\omega) = \int_{E_g}^{\infty} n_s(\hbar\omega) d(\hbar\omega) $$
Dropping the minus one denominator we can write an approximate analytical solution,
$$ N_a \approx f_s\frac{2\pi}{c^2h^3 \beta_s^3}\left(2 + E_g \beta_s \left(2 + E_g \beta_s\right)\right) \exp(-E_g\beta_s) $$
where $\beta_s = 1/\left(kT_s\right)$.
Let's assume our solar cell is Silicon with bandgap $E_g = 1.1\text{eV}$.
Then $N_a \approx 1\times10^{22}$ photons m$^{-2}$s$^{-1}$.
We can use the generalised <PERSON> expression to describe how a solar cell emits light as a function of the applied bias voltage $V$,
$$ n_e(\hbar\omega, V) = \left(\hbar\omega\right)^2\frac{2\pi}{c^2h^3}\frac{1}{\exp{\left( \frac{\hbar\omega - qV}{kT_e} \right)}-1} $$
We want to integrate this to know the total number of photons emitted,
$$ N_e = \int_{0}^{\infty} a(\hbar\omega) n_e(\hbar\omega, V) d(\hbar\omega) = \int_{E_g}^{\infty} n_e(\hbar\omega, V) d(\hbar\omega) $$
again dropping the minus one we can write,
$$ N_e(V) \approx \frac{2\pi}{c^2h^3 \beta_e^3}\left(2 + E_g \beta_e \left(2 + E_g \beta_e\right)\right) \exp(-\left(E_g - qV\right)\beta_e) $$
where the temperature of the solar cell is $T_e=300\text{K}$ and $\beta_e = 1/\left(kT_e\right)$.
We are now in the position to calculate the current voltage characteristics of an ideal solar cell. 246