Multiparticle azimuthal correlations of prompt
D
0
mesons are measured in Pb-Pb collisions at a nucleon-nucleon center-of-mass energy of
√
s
NN
=
5.02
 
 
TeV
. For the first time, a four-particle cumulant method is used to extract the second Fourier coefficient of the azimuthal distribution (
v
2
) of
D
0
mesons as a function of event centrality and the
D
0
transverse momentum. The ratios of the four-particle
v
2
values to previously measured two-particle cumulant results provide direct experimental access to event-by-event fluctuations of charm quark azimuthal anisotropies. These ratios are also found to be comparable to those of inclusive charged particles in the event. However, hints of deviations are seen in the most central and peripheral collisions. To investigate the origin of flow fluctuations in the charm sector, these measurements are compared to a model implementing fluctuations of charm quark energy loss via collisional or radiative processes in the quark-gluon plasma. These models cannot quantitatively describe the data over the full transverse momentum and centrality ranges, although the calculations with collisional energy loss provide a better description of the data.