  The aim of this article is to understand the role of dust reference fields,
often also called clocks, on cosmological perturbations around a flat FLRW
universe. We derive the Mukhanov-Sasaki (MS) equation for the Brown-Kuchar (BK)
and Gaussian dust models, which both consider four dust fields as reference
fields. The reduced phase space of Dirac observables, that is the
gauge-invariant part of the theory, is constructed by means of an observable
map applied to all elementary phase space variables of the coupled system,
consisting of gravity, a massive scalar field and the dust degrees of freedom.
The evolution of these observables is governed by a so called physical
Hamiltonian which can be derived once the set of reference fields are chosen
and differs for each model. First, the reduced phase space as well as the
corresponding equations of motion are derived for full general relativity. Then
from this, the gauge-invariant version of the equations of motion for the
background are derived which contain a fingerprint of the dust reference
fields. Afterwards we study linear cosmological perturbations around a FLRW
metric using the scalar-vector-tensor decomposition and derive the equation of
motion for the MS variable in this formalism for a chosen set of variables on
the reduced phase space and expressed in terms of Dirac observables. The MS
equation involves additional contributions that can be understood as back
reactions from the dust reference fields. These additional dust contributions
to the MS equation were absent if the dust energy and momentum density as well
as their perturbations are vanishing. The nature of the correction terms
suggest that BK and Gaussian dust reference fields contribute differently.
Using numerical simulations we study the behavior of the dust contributions to
the MS equation during inflation.
