PHYSICAL REVIEW LETTERS 128, 134101 (2022)

Classical Physics and Blackbody Radiation

Jiao Wang ,1,2 Giulio Casati,3,4 and Giuliano Benenti 3,5,6
1
Department of Physics and Key Laboratory of Low-Dimensional Condensed Matter Physics
(Department of Education of Fujian Province), Xiamen University, Xiamen 361005, Fujian, China
2
Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou 730000, Gansu, China
3
Center for Nonlinear and Complex Systems, Dipartimento di Scienza e Alta Tecnologia,
Università degli Studi dell’Insubria, via Valleggio 11, 22100 Como, Italy
4
International Institute of Physics, Federal University of Rio Grande do Norte, Campus Universitário—Lagoa Nova,
CP. 1613, Natal, Rio Grande Do Norte 59078-970, Brazil
5
Istituto Nazionale di Fisica Nucleare, Sezione di Milano, via Celoria 16, 20133 Milano, Italy
6
NEST, Istituto Nanoscienze-CNR, I-56126 Pisa, Italy
(Received 15 June 2021; accepted 10 March 2022; published 31 March 2022)

We investigate the properties of the blackbody spectrum by direct numerical solution of the classical
equations of motion of a one-dimensional model that contains the essential general features of the field-
matter interaction. Our results, which do not rely on any statistical assumption, show that the classical
blackbody spectrum exhibits remarkable properties: (i) a quasistationary state characterized by scaling
properties, (ii) consistency with the Stefan-Boltzmann law, and (iii) a high-frequency cutoff. Our Letter is a
preliminary step in the understanding of statistical properties of infinite-dimensional systems.

DOI: 10.1103/PhysRevLett.128.134101

Introduction.—The inability of classical physics to formally in contradiction with the Rayleigh-Jeans law,
account for the experimentally observed frequency spec- according to which E ∝ T.
trum of blackbody radiation is at the origin of quantum What about the solution of the nonlinear classical
theory. In spite of desperate attempts, the falloff of the Newton-Maxwell equations of motion that govern the
blackbody curve at high frequencies could not be explained field-matter interaction? Will dynamics agree, and to what
by classical mechanics. To be more precise, however, this is extent, with the above classical statistical and thermo-
not a failure of classical theory per se. Indeed, it is the dynamics predictions? After 150 years, this question
classical theory with the additional assumption of energy concerning a fundamental problem in the development
equipartition which leads to the Rayleigh-Jeans radiation of modern physics remains unsolved.
law, implying the unphysical ultraviolet catastrophe [1]. Here, we numerically integrate the exact classical equa-
On the other hand, classical ergodic theory, from which tions of motion of a blackbody model, without any
equipartition theorem follows, is valid only for systems statistical assumption, and we show that classical dynamics
with a finite number of degrees of freedom. Ergodicity, that leads to a quasistationary state that is consistent with the
is, independence of time averages on initial conditions, Stefan-Boltzmann law and where the energy distribution
does not imply equipartition of energy for a system with over the field normal modes exhibits an exponentially
infinite degrees of freedom, since there is no invariant decreasing tail. Therefore, the solution of the classical
measure at hand to define a microcanonical ensemble. A equations of motion naturally leads to a high-frequency
main difficulty here stems from the fact that the radiation cutoff of the blackbody spectrum.
field has N m ¼ ∞ degrees of freedom (modes) and that the Model.—To investigate the dynamics of charged par-
two limits N m → ∞ and time t → ∞ do not commute. ticles interacting with the electromagnetic field in a cavity
Therefore, the question of what classical mechanics pre- is a formidable task. Here we consider a variant of a model
dicts on the properties of the radiation field in equilibrium introduced long ago [3] and afterward studied in several
with matter is still open. papers [4–10]. In spite of its simplicity, our model retains
Even more interesting is the Stefan-Boltzmann law, the general essential features of field-matter interaction:
which states that at temperature T the total radiation energy (i) the field modes can only exchange energy via interaction
E ∝ T 4 . It is remarkable that this formula, which is well in with matter and (ii) the free electromagnetic field in the
agreement with experimental data, was derived by cavity is by itself a linear system and nonlinearity is
Boltzmann in 1884 on a purely classical thermodynamics provided by matter.
basis [2]. As such, it should be a consequence of the A sketch of our model is drawn in Fig. 1. It consists of an
classical equations of motion. Notice that this result is electromagnetic field confined in between two parallel,

0031-9007=22=128(13)=134101(5) 134101-1 © 2022 American Physical Society

 PHYSICAL REVIEW LETTERS 128, 134101 (2022)
P
time. In Eq. (1), the term ð1=2mÞðP1 − ε k ak qk Þ2
accounts for the interaction betweenpthe ffiffiffiffiffiffiffi charged plate
and the normal modes, where ε ¼ 2σ π=l is the matter-
field coupling parameter. We assume that fðxÞ is an even
function, and therefore, even normal modes do not interact
with the charged plate, whileR for an odd normal mode with
wave number 2k − 1, ak ¼ −δ δ fðxÞ cosðω x=cÞdx, where
k
c is the speed of light and ωk ¼ ðπc=2lÞð2k − 1Þ.
If the charge is removed, ε ¼ 0, the plates and the modes
are decoupled. We choose, as interaction potential among
the matter degrees of freedom, the lattice ϕ4 model [11,12]

1 1
Ṽðzj Þ ¼ γz4j ; Vðzj−1 ; zj Þ ¼ κðzj − zj−1 Þ2 : ð2Þ
4 2

We checked that, for N p > 4 and for an average energy per

FIG. 1. Schematic drawing (top) of the classical radiant cavity plate larger than 0.1 (in our units l ¼ π, c ¼ m ¼ σ ¼ 1,
model. A charged plate of width 2δ is placed at the center of the and the Boltzmann constant kB ¼ 1; moreover, we set
cavity bounded by two mirrors a distance 2l apart. The charged
γ ¼ κ ¼ 1), this model is chaotic with N p positive
plate interacts with a fixed number of neutral plates (to the left in
the figure) and with the normal modes of the electromagnetic Lyapunov exponents.
field. The topology of interactions is also shown (bottom), where Initially, all the energy is assigned to matter only (the
the gray dot, red dots, and blue dots are for the charged plate, the plates), and we focus on how the energy is transferred and
neutral plates, and the field modes, respectively. distributed among the field normal modes. For the charge
distribution we choose fðxÞ ¼ A expðδ2 =ðx2 − δ2 ÞÞ (jxj < δ)
(a compactly supported C∞ function, with A normalization
perfectly reflecting plane mirrors, a distance 2l apart. We
constant, A ¼ 45.04… for the width parameter δ ¼ 0.05
take Cartesian coordinates xyz with the x axis normal to the
used in our simulations). We have verified that, qualita-
mirrors and restrict to excitations only dependent on x, thus
tively, the results do not depend on the particular choice of
getting a one-dimensional radiant cavity, the normal modes
the charge distribution, provided it is a smooth function
of which have angular frequencies ωk . We then introduce
[13]. In addition, we take N p ¼ 16 and we have verified
N p coupled plates, which are all parallel to the mirrors and
that the whole system remains chaotic during the entire
which are constrained to move only in the z direction. Only
relaxation process accessible to simulations (t ∼ 109 ) (see
the plate positioned midway between the mirrors is charged
Supplemental Material [13]).
and therefore interacts with the normal modes of the field.
Results.—Since our simulations are for energies well
We denote by σ and m the charge and mass densities per
above the stochasticity threshold, for any finite number of
unit surface of the plate, and fðxÞ is the normalized
field modes, equipartition is expected among the degrees of
(transverse) distribution of charge in the plate, whose
freedom. Our first step is therefore to investigate how long
thickness is 2δ.
it takes for equipartition to set in. As shown in Fig. 2,
The Hamiltonian of the full system, plates plus field, can
equipartition is reached among the plates after a relatively
be written as
short time τp that, basically, does not depend on the number
Np  2  N m of field modes. Instead, global equipartition among
X Pj
H¼ þ Ṽðzj Þ þ Vðzj−1 ; zj Þ þ Ṽðz1 Þ field and matter degrees of freedom is reached after a time
j¼2
2m τm that exponentially increases with N m . The best fitting
 X 2 suggests τm ∼ e1.29N m , in clear contrast with τp ≈ 103 [13].
1 ∞
1X ∞
þ P1 − ε ak qk þ ðp2 þ ω2k q2k Þ; ð1Þ If we extrapolate this exponential dependence of the
2m k¼1
2 k¼1 k equipartition time τm to larger N m values, then it turns
out that already for N m ≈ 48 field normal modes, a time of
where Pj and zj are the conjugate momentum and the order of the age of the Universe is required in order to
displacement of the jth plate, Ṽ and V are the on site reach equipartition (in a cavity of l ∼ 1 m). Note that the
and interacting potentials of plates, and pk and qk are two limits t → ∞ and N m → ∞ do not commute. If one
conjugate variables of the kth normal mode of the field with takes the limit t → ∞ first, then one will always get
frequency ωk. Note that the field normal modes are infinite equipartition, while for the blackbody problem it is
but, as we shall see below, only a finite, small number necessary to take the limit N m → ∞ first. In order to
N m ðtÞ of them are actually involved during the computation simulate this latter situation, in our computations we take

134101-2
 PHYSICAL REVIEW LETTERS 128, 134101 (2022)

hEm;k ðtÞi of the normal modes over their frequency. The

main feature of the distribution hEm;k i is the presence of a
plateau followed by a rapid, exponential decay. The
average value hEm;k i of the plateau normal mode energies
is equal to the average energy hEp;j i of the plates. This
suggests that the field modes in the plateau are thermalized
with matter, at temperature T ¼ hEp;j i. This thermalization
process evolves in time very slowly; indeed, as seen in
Fig. 3(a), it takes about 3 orders of magnitude in time (from
105 to 108 ) in order to increase only by four the number of
thermalized modes. This observation is consistent with the
results of Fig. 2 and suggests that the system relaxes to a
quasistationary state, after which it evolves very slowly,
logarithmically in time.
FIG. 2. The relaxation time of plates (τp ) and of both plates and
field modes (τm ), where N m field modes are considered in the
It is remarkable that the quasistationary state has a clear
simulation. Here the number of plates N p ¼ 16, and the total scaling property. This is shown in Figs. 3(b) and 3(c),
energy Etot ¼ 24 is initially assigned to the charged plate as where we plot hEm;k ðtÞi=TðtÞ over ½ωk − α lnðtÞ and over
kinetic energy. For other initial conditions where the total energy β½ωk − α lnðtÞ, where α and β depend on the total energy
is distributed randomly among the plates, τm ∼ e1.29N m still holds. only, α ¼ 0.23E0.24tot and β ¼ 3.22Etot
−0.22
based on the best
fitting of the data for 10 ≤ t ≤ 10 and 24 ≤ Etot ≤ 210 .
4 8

We cannot give a rigorous explanation for the appear-

N m ¼ 120, which, in view of the results of Fig. 2, is
ance of an exponential cutoff. Nevertheless, we can provide
sufficiently large to neglect the higher normal modes.
an intuitive understanding by observing (see Fig. 4) that, as
We now turn to the energy distribution over the field
is typically the case for nonlinear interacting systems, the
normal modes. In Fig. 3(a) we plot the average energy
power spectrum of the charged plate has a finite frequency
band, followed by an exponential tail. For the field modes
within this band, thermalization with matter is achieved
(a)
rapidly. On the other hand, for the field modes with
frequencies in the exponential tail, thermalization requires
exponentially long timescales. We can see from Fig. 4 that
the power spectrum is, in practice, independent of the
number N p of plates (for a fixed energy per plate), while for
a given N p the bandwidth increases with the total energy
Etot . Such dependence is consistent with the reduction of
the thermalization times (up to a given frequency) with
increasing Etot .
(b) (c)

FIG. 3. (a) The average energy of the plates (full symbols) and
of the field modes (half filled symbols), for total energy Etot ¼ 24 (a) (b)
at different times. (b) hEm;k i scaled by the temperature of plates
and the mode frequencies being shifted by αln(t). (c) hEm;k i
versus the scaled and shifted frequencies at a given time. The
results for different total energies collapse on the same curve. The FIG. 4. Power spectrum of the time series P1 ðtÞ for the
energy of the system is initially randomly assigned to the plates as momentum of the charged plate, in a system of N p plates with
kinetic energy. total energy Etot, for (a) Etot =N p ¼ 8 and (b) N p ¼ 16.

134101-3
 PHYSICAL REVIEW LETTERS 128, 134101 (2022)

FIG. 5. The pressure p, the energy density u of the cavity, and

the temperature T of the plates as functions of time.

FIG. 6. The energy density u in the cavity versus the temper-

Note that the total energy of the modes in the exponential ature T of the plates at various evolution times. At a given
tail is negligible compared with that of the thermalized time, the seven data points correspond to total energies
modes on the plateau. Therefore, it is tempting to con- Etot ¼ 24 ; 25 ; …; 210 , respectively. At each time, the data exhibit
jecture that the quasistationary state should be approxi- a power law relation u ∝ T B , and the best fitting shows that the
mately described by equilibrium thermodynamics. To this power law exponent BðtÞ slowly changes from 1.07 to 1.21 as
end, following Boltzmann, we start from the fundamental time changes by 4 orders of magnitude. The fitting line in the
thermodynamic relation dU ¼ TdS − pdV, which implies inset gives the expression 2 − BðtÞ ∼ ðln tÞ−0.33 .
 
∂U  ∂p inset of Fig. 6 suggests that the value 2 might be
¼ −p þ T  ; ð3Þ
∂V T ∂T V approached logarithmically in time.
Summary and discussion.—In summary, we have studied
where U ¼ UðT; VÞ is the total energy, S is the entropy, p a model of a classical radiant cavity. If one considers a
is the pressure of the radiation field, V is the volume of the fixed, finite number N m of field modes, then the system, as
cavity, and T is the temperature of the plates. The relation expected, approaches equipartition, even though the relax-
between the energy of the electromagnetic radiation and its ation times increase exponentially with the number N m of
pressure has been derived by Boltzmann and for the one- modes. If instead one takes the limit N m → ∞ first, as it is
dimensional case reads required in order to study the blackbody radiation, then the
system relaxes to a quasistationary state characterized by a
p ¼ u; ð4Þ low frequency plateau, followed by a high frequency
exponential cutoff. The computed field energy density
where u ¼ U=V is a universal function of temperature and turns out to be consistent with the Stefan-Boltzmann law.
independent of volume. From Eq. (3) one then has Our model contains the basic ingredients of an electro-
magnetic field interacting with matter. Since interaction
uðTÞ ¼ CT 2 ; ð5Þ among the field modes is mediated by mechanical degrees
of freedom, in general, we expect thermalization of the field
where C is a constant. This is the Stefan-Boltzmann law in modes with matter to be effective only within the frequency
one dimension. Quite remarkably, these predictions based bandwidth of mechanical motion, with slow thermalization
on equilibrium thermodynamics are consistent with our outside such bandwidth. We therefore conjecture that the
model in the quasistationary state. In Fig. 5, we plot p, u, appearance of a time-dependent cutoff is a general feature of
and T versus the evolution time for Etot ¼ 24. Here p is classical dynamics of matter-field interaction. Nevertheless, it
numerically computed (see Supplemental Material [13]) would be interesting to investigate higher-dimensional mod-
and T ¼ hEp;j i is the temperature of matter. It can be els or, more generally, other classical field theories. Such
noticed that the relation (4) between energy density and studies could help our understanding of the statistical proper-
pressure sets in almost immediately. In Fig. 6, we plot u ties of nonlinear dynamical systems with infinite degrees of
versus T for various values of the total energy of the freedom. In this frame, the results presented here on a
system. It is seen that u ∝ T B where the exponent BðtÞ classical model can be considered a preliminary step before
increases from 1.07 to 1.21, as time t goes from 104 to 108 . addressing ergodicity in quantum field theories. This problem
This dependence is compatible with Eq. (5). Indeed the will require a nontrivial extension of concepts and tools

134101-4
 PHYSICAL REVIEW LETTERS 128, 134101 (2022)

recently developed for the investigation of thermalization and [8] R. Livi, M. Pettini, S. Ruffo, and A. Vulpiani, J. Phys. A 20,
localization in many-body quantum systems [14–18]. 577 (1987).
[9] C. Alabiso, M. Casartelli, and S. Sello, J. Stat. Phys. 54, 361
We acknowledge support by the NSFC (Grants (1989).
No. 12075198 and No. 12047501) and by the INFN [10] G. Benenti, G. Casati, and I. Guarneri, Europhys. Lett. 46,
through the project QUANTUM. 307 (1999).
[11] K. Aoki and D. Kusnezov, Phys. Lett. A 265, 250
(2000).
[12] G. Parisi, Europhys. Lett. 40, 357 (1997).
[1] F. R. S. Lord Rayleigh, Philos. Mag. 49, 539 (1900); J. H. [13] See Supplemental Material at http://link.aps.org/supplemental/
Jeans, The Dynamical Theory of Gases, 4th ed. (Cambridge 10.1103/PhysRevLett.128.134101, where Ref. [10] is in-
University Press, Cambridge, England, 2009), 10.1017/ cluded, for results obtained for other charge distributions
CBO9780511694370. and details on numerical simulations.
[2] L. Boltzmann, Ann. Phys. (Berlin) 258, 291 (1884); Nature [14] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore,
(London) 51, 413 (1895). Rev. Mod. Phys. 83, 863 (2011).
[3] P. Bocchieri, A. Crotti, and A. Loinger, Lett. Nuovo [15] R. Nandkishore and D. A. Huse, Annu. Rev. Condens.
Cimento 4, 741 (1972). Matter Phys. 6, 15 (2015).
[4] P. Bocchieri, A. Loinger, and F. Valz-Gris, Nuovo Cimento [16] F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G.
B 19, 1 (1974). Zelevinsky, Phys. Rep. 626, 1 (2016).
[5] G. Casati, I. Guarneri, and F. Valz-Gris, Phys. Rev. A 16, [17] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv.
1237 (1977). Phys. 65, 239 (2016).
[6] G. Benettin and L. Galgani, J. Stat. Phys. 27, 153 (1982). [18] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.
[7] G. Casati, I. Guarneri, and F. Valz-Gris, J. Stat. Phys. 30, Mod. Phys. 91, 021001 (2019).
195 (1983).

134101-5