  A recent reformulation [1] of the problem of Coulomb gases in the presence of
a dynamical dielectric medium showed that finite temperature simulations of
such systems can be accomplished on the basis of completely local Hamiltonians
on a spatial lattice by including additional bosonic fields. For large systems,
the Monte Carlo algorithm proposed in Ref. [1] becomes inefficient due to a low
acceptance rate for particle moves in a fixed background multiboson field. We
show here how this problem can be circumvented by use of a coupled
particle-multiboson update procedure that improves acceptance rates on large
lattices by orders of magnitude. The method is tested on a one-component plasma
with neutral dielectric particles for a variety of system sizes.
