Magnetars are neutron stars (NSs) with extreme magnetic fields1 of strength 5 × 1013−1015 G. These fields are generated by dynamo action during the proto-NS phase, and are expected to have both poloidal and toroidal components2–6, although the energy of the toroidal component could be ten times larger7. Only the poloidal dipolar field can be measured directly, via NS spin-down8. The magnetic field provides heating and governs how this heat flows through the crust9. Magnetar thermal X-ray emission in quiescence is modulated with the rotational period of the NS, with a typical pulsed fraction 10–58%, implying that the surface temperature is substantially non-uniform despite the high thermal conductivity of the star’s crust. Poloidal dipolar fields cannot explain this large pulsed fraction10,11. Previous two-dimensional simulations12,13 have shown that a strong, large-scale toroidal magnetic field pushes a hot region into one hemisphere and increases the pulsed fraction. Here, we report three-dimensional magneto-thermal simulations of magnetars with strong, large-scale toroidal magnetic fields. These models, combined with ray propagation in curved spacetime, accurately describe the observed light curves of 10 out of 19 magnetars in quiescence and allow us to further constrain their rotational orientation. We find that the presence of a strong toroidal magnetic field is enough to explain the strong modulation of thermal X-ray emission in quiescence.