  We study the leading-order nonlinear gravitational waves (GWs) produced by an
electromagnetic (EM) stress in reheating magnetogenesis scenarios. Both
nonhelical and helical magnetic fields are considered. By numerically solving
the linear and leading-order nonlinear GW equations, we find that the GW energy
from the latter is usually larger. We compare their differences in terms of the
GW spectrum and parameterize the GW energy difference due to the nonlinear
term, $\Delta\mathcal{E}_{\rm GW}$, in terms of EM energy $\mathcal{E}_{\rm
EM}$ as $\Delta\mathcal{E}_{\rm GW}=(\tilde p\mathcal{E}_{\rm EM}/k_*)^3$,
where $k_*$ is the characteristic wave number, $\tilde p=0.84$ and $0.88$ are
found in the nonhelical and helical cases, respectively, with reheating around
the QCD energy scale, while $\tilde p=0.45$ is found at the electroweak energy
scale. We also compare the polarization spectrum of the linear and nonlinear
cases and find that adding the nonlinear term usually yields a decrease in the
polarization that is proportional to the EM energy density. We parameterize the
fractional polarization suppression as $|\Delta \mathcal{P}_{\rm
GW}/\mathcal{P}_{\rm GW}|=\tilde r \mathcal{E}_{\rm EM}/k_*$ and find $\tilde r
= 1.2 \times 10^{-1}$, $7.2 \times 10^{-4}$, and $3.2 \times 10^{-2}$ for the
helical cases with reheating temperatures $T_{\rm r} = 300 {\rm TeV}$, $8 {\rm
GeV}$, and $120 {\rm MeV}$, respectively. Prospects of observation by pulsar
timing arrays, space-based interferometers, and other novel detection proposals
are also discussed.
