  We study a model where two scalar fields, that are subdominant during
inflation, decay into radiation some time after inflation has ended but before
primordial nucleosynthesis. Perturbations of these two curvaton fields can be
responsible for the primordial curvature perturbation. We write down the full
non-linear equations that relate the primordial perturbation to the curvaton
perturbations on large scales, and solve them in a sudden-decay approximation.
We calculate the power spectrum of the primordial perturbation, and finally go
to second order to find the non-linearity parameter, fNL.
  Not surprisingly, we find large positive values of fNL if the energy
densities of the curvatons are sub-dominant when they decay, as in the single
curvaton case. But we also find a novel effect, which can be present only in
multi-curvaton models: fNL becomes large even if the curvatons dominate the
total energy density in the case when the inhomogeneous radiation produced by
the first curvaton decay is diluted by the decay of a second nearly homogeneous
curvaton. The minimum value min(fNL)=-5/4 which we find is the same as in the
single-curvaton case.
  Using (non-)Gaussianity observations, Planck can be able to distinguish
between single-field inflation and curvaton model. Hence it is important to
derive theoretical predictions for curvaton model. From particle physics point
of view it is more natural to assume multiple scalar fields (rather than just
one ``curvaton'' in addition to inflaton). Our work updates the theoretical
predictions of curvaton model to this case.
