\section*{Real-time polarimetry of hyperpolarized \textsuperscript{13}C nuclear spins using an atomic magnetometer}
\subsection*{Supporting Information}

\vspace{0.1cm}

\vspace{0.3cm}

\noindent \textbf{Kostas Mouloudakis\textsuperscript{a}, Marc Azagra\textsuperscript{b}, S  Bodenstedt\textsuperscript{a}, Irene Marco-Rius\textsuperscript{b}, Morgan W. Mitchell\textsuperscript{a,c}, Michael C. D. Tayler\textsuperscript{a,*}}\\

\noindent \emph{\textsuperscript{a} \quad \footnotesize ICFO-Institut de Ci\`encies Fot\`oniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain.}

\noindent \emph{\textsuperscript{b} \quad \footnotesize IBEC -- Institute for Bioengineering of Catalonia, 08028 Barcelona, Spain}

\noindent \emph{\textsuperscript{c} \quad \footnotesize ICREA -- Instituci\'{o} Catalana de Recerca i Estudis Avan\c{c}ats, 08010 Barcelona, Spain.}

\noindent \emph{\textsuperscript{*} \quad \footnotesize Corresponding author.  Email: michael.tayler@icfo.eu}

\vspace{0.2cm}

\subsection*{Contents}
\nameref{sec:setup} \dotfill \pageref{sec:setup}\\
\nameref{sec:shield} \dotfill \pageref{sec:shield}\\
\nameref{sec:opm} \dotfill \pageref{sec:opm}\\
\nameref{sec:solenoid} \dotfill \pageref{sec:solenoid}\\
\nameref{sec:pulses} \dotfill \pageref{sec:pulses}\\
\nameref{sec:misc} \dotfill \pageref{sec:misc}\\
\newpage

\subsection*{S.1 \,\, Experimental detail} \label{sec:setup}

The magnetic field of the nuclear-spin-polarized liquid sample was detected using the experimental setup illustrated in \autoref{fig:Setup_photographs} and \autoref{fig:Setup_schematic}.  Main parts of the setup are described in the following subsections.
\begin{figure}[h]%
    \centering
	\includegraphics[width=.7\textwidth]{photos.png}
	\caption{Photographs of the experimental apparatus in use at the molecular-imaging-for-precision-medicine group at the Catalonia Institute for Bioengineering (IBEC): (a) the nuclear spin polarimeter, on top of a table approximately 1 m away from a commercial dissolution-DNP polarizer.  The main visible feature of the setup is the cylindrical MuMetal magnetic shield; (b) a view of the interior of the magnetic shield, containing the coils for XY4 pulses on a flexible printed circuit board, the magnetometer, and a solenoid coil that surrounds the sample tube containing hyperpolarized liquid.  A schematic version of (b) appears in the lower part of \autoref{fig:Setup_schematic} .}
	\label{fig:Setup_photographs}
\end{figure}
\begin{figure}[h]%
    \centering
	\includegraphics[width=.7\textwidth]{setup-schematic.png}
	\caption{Schematic view of the experimental setup for real-time nuclear spin polarimetry.  The top half of the figure illustrates the magnetically shielded portion of the setup from a vertical point of view, down the $z$ axis. The bottom half shows the view along the $y$ axis, parallel to the laser beam.}
	\label{fig:Setup_schematic}
\end{figure}

\subsection*{S.1.1 Magnetic shielding} \label{sec:shield}

A commercial magnetic shield was used to passively screen background magnetic fields in the laboratory, up to several mT in magnitude, down to a residual of between 10 and 30 nT.  The shield comprised three concentric outer layers of MuMetal and one innermost layer of ferrite (MS-1LF, Twinleaf LLC).  The shield as receiver from the manufacturer also contained a set of XYZ printed-circuit coils on a flexible polyimide film attached to the interior of the ferrite layer, as shown in the \autoref{fig:Setup_photographs}(b).  These coils are normally intended to be used to supply compensatory magnetic fields along $x$-, $y$- and $z$-axes to achieve zero residual magnetic field inside the shield.  In the present work, the $x$- and $y$-axis coils were used to supply the XY4 pulses (indicated in Figure 1 of the main manuscript) and the $z$-axis coil (Z) was unused.  Field-to-current ratios were 69 nT/mA for the X coil and 130 nT/mA for the Y coil.

\subsection*{S.1.2 Atomic magnetometer}\label{sec:opm}

A linearly polarized 1 mW laser beam was generated by a vertical-cavity surface-emitting laser at the \textsuperscript{87}Rb $\rm{D}_1$ transition wavelength, around 795 nm.  The beam was passed through a Keplerian telescope to adjust its width to around 3 mm and then a zero-order quarter-wave plate to generate circularly polarized light.  The light then passed into the magnetic shield through the end hole and then through an optically dense alkali medium enclosed in a MEMS vapor cell.  Cell external dimensions were 5 mm by 10 mm by 1.5 mm and internal path length $1$ \si{\milli\meter}.  The cell was locally heated to a temperature of around \SI{150}{\celsius} using a resistive heater.  The beam passing through the vapor was detected outside the magnetic shield, where it was focused on a low-noise amplified photodetector (Thorlabs model PDA36A2) and the resulting output voltage sampled using a 24-bit analog-to-digital converter integrated circuit at 48 ksps (Cirrus Logic CS4272) interfaced with a microcontroller (ARM Cortex M4F).

To operate the above apparatus as a magnetometer, the local magnetic field at the atomic vapor cell was modulated by applying a sinusoidal field along the $z$ direction of 20 nT in amplitude and approximately 500 Hz in frequency. The field was supplied using a miniature biplanar printed-circuit coil set\cite{Tayler2022PRApplied} positioned around the cell oven, as shown in \autoref{fig:Setup_photographs}(b) in green.  During modulation, the photodetector signal was sampled using the CS4272 and demodulated at the same frequency, then digital lowpass filtered to produce a quadrature signal whose amplitude was proportional to the dc bias field $B_z$.  This field was the sum of the residual field inside the magnetic shield and the nuclear magnetic field of the hyperpolarized sample, seen at the position the atoms.

\subsection*{S.1.3 Piercing solenoid} \label{sec:solenoid}

A 10-inch-long, 10-mm-o.d., 8.5-mm-i.d. flat-bottomed glass tube was used to contain the NMR sample during the magnetometry measurements.  This was wrapped along its length with 0.15 mm diameter enameled copper wire to form a single layer solenoid coil, which when energized supplied a \SI{30}{\micro\tesla} bias field along the $z$ axis; this field ensured that the spin polarization was always parallel to $z$, even when the sample was transported in and out of the magnetic shield.  The bottom of the tube/solenoid structure sat atop the magnetometer oven as shown in \autoref{fig:Setup_photographs}(b).

\subsection*{S.1.4 Pulsed field control} \label{sec:pulses}

The alternating $(\pi)_x \equiv (\pi/2)_x(\pi)_y(\pi/2)_x$ and $(\pi)_y \equiv (\pi/2)_y(\pi)_{-x}(\pi/2)_y$ composite pulses constituting the XY4 pulse sequence shown in Figure 1 of the main paper were applied via bidirectional switching of direct current in the X and Y coils, as shown schematically in \autoref{fig:Setup_schematic}.  Current switching in both of the coils was performed using a dual-channel H-bridge circuit (Toshiba TB6612FNG on Pololu carrier board).  The logic lines of the H-bridge circuit were controlled using the M4F microcontroller's digital input-output and were asynchronous with CS4272 data acquisition.  The approximate duration of each $\pi$ composite pulse was \SI{50}{\micro\second} on \textsuperscript{1}H for a peak current of 0.5 A.

\subsection*{S.1.5 Additional notes} \label{sec:misc}

The table-top setup involving the magnetic shield, magnetometer and control electronics was constructed at ICFO and transported to IBEC in the rear passenger seat of a car belonging to one of the authors.  The experiments described in the main manuscript were then performed at IBEC in the same laboratory as a commercial dissolution-DNP polarizer (HyperSense, Oxford Instruments Ltd.), as shown in \autoref{fig:Setup_photographs}(a), on a vibrationally undamped table located approximately 1 m away from the polarizer magnet.  The distance was chosen so that the magnetometer would be as close as possible to the polarizer without impeding access to its working parts.

\section*{KEYWORDS}
\noindent
Nuclear magnetism;
Hyperpolarization;
Optical Magnetometry;
Polarimetry;
Quantum control.

\section*{Abstract}
\noindent We introduce a method for non-destructive quantification of nuclear spin polarization, of relevance to hyperpolarized spin tracers widely used in magnetic resonance from spectroscopy to in vivo imaging.  In a bias field of around \SI{30}{\nano\tesla} we use a high-sensitivity miniaturized \textsuperscript{87}Rb vapor magnetometer to measure the field generated by the sample, as it is driven by a windowed dynamical decoupling pulse sequence that both maximizes the nuclear spin lifetime and modulates the polarization for easy detection. We demonstrate the procedure applied to a \SI{0.08}{\molar} hyperpolarized [1--\textsuperscript{13}C]-pyruvate solution produced by dissolution dynamic nuclear polarization, measuring polarization repeatedly during natural decay at Earth's field.  Application to real-time quality monitoring of hyperpolarized substances is discussed.

\section*{Main text}

\noindent Magnetic resonance spectroscopy (NMR) and imaging (MRI) of nearly fully polarized (also known as ``hyperpolarized'') nuclear-spin ensembles is an emerging methodology\cite{Aquaro2014,Comment2014Biochem,Brindle2015,Wang2019,Jrgensen2022} for probing cellular biochemistry
in vitro\cite{},
in vivo\cite{MarcoRius2017emagres,MarcoRius2020MRMPBM,vonMorze2017MRM,Can2022CommBio
and in perfused organs \cite{Merritt2007PNAS}.
The magnetic moment of a hyperpolarized substance, which is proportional to the nuclear polarization, can be three to four orders of magnitude larger than that of thermally polarized substances\cite{ArdenkjrLarsen2003PNAS}, resulting in MRI pixels that contain chemically specific signals high above the noise floor even at low concentration.
In pre-clinical\cite{Ros2021STAR} and clinical\cite{Gallagher2020PNAS,Woitek2021BJC} monitoring of metabolic diseases, hyperpolarized small molecules such as pyruvic acid are typically injected into a biosystem where metabolism takes place.
Dynamic MRI of the downstream metabolites then provides biodistribution maps and kinetic rates to distinguish between healthy and diseased tissue.

In these and other applications, such as noble-gas MRI, there are several motives to quantify nuclear spin hyperpolarization.  The magnitude or degree of polarization, $P$ (a dimensionless number between zero and one), determines the effective concentration of the substrate seen by MRI and therefore is a crucial parameter in kinetic models of spin-label transfer along a metabolic pathway, especially in cases where hyperpolarized agents are administered to the subject repeatedly\cite{Lees2020NMRBiomedicine}.  Additionally, any exogenous substance introduced in vivo is subject to strict quality control and dose requirements.  Given that many hyperpolarization technologies are at an early development stage and not always consistent, $P$ is currently a control variable; for instance, in dissolution dynamic nuclear polarization,
polarization is monitored both before\cite{Elliott2021,Koptyug2022PCCP24} and after\cite{ArdenkjaerLarsen2011NMRBiomed24,Kurhanewicz2019Neoplasia} sample dissolution.
A reliable measurement of $P$ may also provide a cost-per-unit signal enhancement metric to compare different hyperpolarization methods, which may be influential in future large-scale clinical or other applicability of hyperpolarized magnetic resonance.

A standard polarimetry method is \textit{ex post} comparison of the hyperpolarized-NMR signal amplitude against the NMR signal of the same system after polarization has decayed back to a thermal equilibrium value\cite{ArdenkjaerLarsen2011NMRBiomed24,Kurhanewicz2019Neoplasia}.  This approach can be performed on a variety of NMR spectrometer systems in background fields from mT to T, but is slow and requires several calibrations such as tuning of radiofrequency coils, pulse flip angles, and magnetic field; conventionally, although lower-field instrumentation may be less expensive, calibrations become much lengthier due to reduced-efficiency inductive pickup of the spin precession signals\cite{Korchak2012AMR44}.  In special molecules one can also determine polarization from intensity patterns in spin-coupling multiplets, namely where the hyperpolarized substrate contains a second, ancillary spin label to be used as a ``spy'' nucleus, e.g., [1,2-\textsuperscript{13}C\textsubscript{2}]-pyruvate\cite{MarcoRius2013,IMRPhDthesis,Vuichoud2015JMR260} or [2-\textsuperscript{13}C]-acetate\cite{Elliott2021}.
A third method is through detection of the bulk magnetic field, $\mathbf{B}$, produced by the hyperpolarized spins.  For an ensemble of uncoupled spin-1/2 nuclei,
\begin{equation}
  |\mathbf{B}| = \hbar\mu_0 P C \gamma k\,,
  \label{eq:Mz}
\end{equation}
assuming that $\mathbf{B}$ is parallel to the polarization vector, where $\hbar$ is the reduced Planck constant, $\mu_0$ is the vacuum permeability, $C$ is the number density of spins, $\gamma$ is the gyromagnetic ratio and $k$ is a constant of proportionality that depends on the geometry of the sample and the position where $\mathbf{B}$ is measured.
In the case of a uniformly polarized spherical sample of radius $R$, $k=(R/r)^3/3$ for position $r>R$ along the $\mathbf{B}$ direction.

Most bulk-protonated liquids have $|\mathbf{B}|$ on the order of \si{\pico\tesla} to \si{\nano\tesla} for $R/r \approx 1$, at thermal equilibrium in moderate fields, e.g., 1 T.
A dilute hyperpolarized liquid ($PC \approx P_{\rm thermal}C_{\rm bulk}$) could produce fields of comparable magnitude, while a hyperpolarized neat liquid ($C \approx C_{\rm bulk}$, $P\gg P_{\rm thermal}$) will produce a considerably stronger field.
One way to measure $|\mathbf{B}|$ is through Zeeman frequency shifts of a second, thermally polarized spin species in the liquid\cite{Eichhorn2022jacs}. This, in addition to the above methods, however, relies on high-resolution NMR instrumentation (with which comes high cost and limited mobility) and the destruction of a significant part if not all of the hyperpolarized spin order.  In the following, as an alternative, we demonstrate a direct, nondestructive measurement of $|\mathbf{B}|$ using an optically pumped magnetometer (OPM)\cite{Savukov2005PRL,Tayler2017RSI,Yashchuk2004PRL93}.  To obtain $P$, the factor $k$ in \autoref{eq:Mz} is eliminated using a reference sample of known polarization ($P_{\rm ref}$) and concentration ($C_{\rm ref}$):
\begin{equation}
  P = \left( \frac{\gamma_{\rm ref}}{\gamma} \right) \left( \frac{C_{\rm ref}}{C} \right) \frac{|\mathbf{B}|}{|\mathbf{B}_{\rm ref}|}P_{\rm ref}\,.
  \label{eq:PA}
\end{equation}
Here the reference sample is \textsuperscript{1}H spins in thermally polarized water at \SI{2}{\tesla} (\textsuperscript{1}H\textsubscript{2}O; $C_{\rm ref} = $ 110 \si{\molar}, $P_{\rm ref} \approx 6.8 \times 10^{-6}$).

A suitable pulse sequence to detect $\mathbf{B}$, assumed parallel to the $z$ axis, is shown in \autoref{fig:figure1}.  Uniformly spaced $\pi$ spin flips toggle the sign of the $z$ component, producing a square wave that is distinct from background fields\cite{Ganssle2014ANIE}.
Flips are applied as $\pi$ rotations alternately about $x$ and $y$ laboratory-frame axes making an XY4 cycle with repeating element ($\tau/2-\pi_{\rm X}-\tau-\pi_{\rm Y}-\tau-\pi_{\rm X}-\tau-\pi_{\rm Y}-\tau/2$).  The useful feature of XY4 is the dynamical decoupling of the background field from the spin system at the first-order average Hamiltonian level \cite{Bodenstedt2021jpclett}, which allows us to neglect bias fields with strengths $|\mathbf{B}_\mathrm{bias}| \ll 2\pi/\gamma\tau$.

\begin{figure}%
\includegraphics[width=0.9\columnwidth]{./fig1.png}%
\caption{Experimental sequence.  Upon field switching to near-zero fields, alternating $(\pi)_x \equiv (\pi/2)_x(\pi)_y(\pi/2)_x$ and $(\pi)_y \equiv (\pi/2)_y(\pi)_{-x}(\pi/2)_y$ dc pulses (making the XY4 cycle) periodically invert the nuclear magnetization.  A magnetometer measures the $z$ field component throughout.
}
\label{fig:figure1}\vspace{0.5cm}
\includegraphics[width=\columnwidth]{./H2Odata.png}%
\caption{Measured $z$-axis component of magnetic field ($B_z$, relative to background) during the ultralow-field portion of the sequence shown in \autoref{fig:figure1}.  Single-shot signals are measured for \SI{1}{\centi\meter\cubed} of H\textsubscript{2}O using a \textsuperscript{87}Rb-vapor OPM under the following conditions: (top) XY4 decoupling, sample prepolarized at \SI{2}{\tesla}, (middle) XY4 decoupling, prepolarized at Earth's field, (bottom) free precession, i.e., without XY4 pulses, prepolarized at \SI{2}{\tesla}.  The red curves indicate fitting to square- or sine-modulated monoexponential decay curves with the time constants (top) $T_1 = \SI{3.5}{\second}$ and (bottom) $T_2^\ast = \SI{2.0}{\second}$.

\label{fig:figure2}
\end{figure}

The top panel of \autoref{fig:figure2} shows the NMR signal recorded during XY4 ($\tau=50$ ms, 0.2 ms $\pi$-pulse length) for a 1 \si{\centi\meter\cubed} cylindrical volume of 2-T-prepolarized water, using an alkali-metal-vapor OPM inside a magnetic shield (MS-1LF, Twinleaf LLC).  The OPM measures the $z$-axis component, $B_z$, of total field using quadrature demodulation of the Hanle absorption resonance in a \textsuperscript{87}Rb vapor\cite{Shah2007}.  The vapor is confined to a microelectromechanical system  (MEMS) cell of approximate dimensions 4 mm $\times$ 4 mm $\times$ 1.5 mm at a standoff distance of 3 mm from the sample chamber.  The magnetic shield provides a simple and convenient way to isolate the detector from unknown and/or time-dependent background magnetic fields.  For further technical details of the OPM we refer the reader to the Supporting Information, plus past work\cite{Tayler2022PRApplied}.  The signal shown is processed minimally, where data points during each $\pi$ pulse and up to 5 ms afterwards are excluded (where the magnetometer output saturates) and mains-hum noise is suppressed by applying a 45-to-55-Hz band-stop filter.  The solid red curve: $B_z(t) = (190\,\,\rm pT)\times{\rm sgn}(\sin 20 \pi t)\exp(-t/T_1)$ fits to the remaining data points.  The decay time constant, $T_1=3.5$ s, is that of longitudinal relaxation under dynamical decoupling in the effective zero field, inclusive of the contribution due to imperfections in the $\pi$ pulses.

The middle panel of \autoref{fig:figure2} shows the OPM signal when the starting magnetic field is much lower, \SI{50}{\micro\tesla} or around Earth's field, where there is no major alternation seen above the noise background. As a comparison, the lower panel of \autoref{fig:figure2} shows the signal for 2-T-prepolarized water in the absence of XY4 decoupling.  Here the free-decay signal is an exponentially decaying cosine wave $B_z(t) = (190\,\,\rm pT)\times\cos (2.48 \pi t)\exp(-t/T_2^\ast)$ with $T_2^\ast = \SI{2.0}{\second}$, corresponding to precession in the approximately \SI{30}{\nano\tesla} remanent field of the magnetic shield.  The precession signal undergoes faster decoherence ($T_2^\ast<T_1$) due to gradients in the precession field, which are otherwise decoupled under XY4.

\begin{figure*}%
\includegraphics[width=0.95\textwidth]{./fig3new.png}%
\caption{The NMR signals recorded for \SI{0.08}{\molar} hyperpolarized [1-\textsuperscript{13}{C}]-pyruvate under \textsuperscript{13}{C} XY4, $\tau=50$ ms. Upper panel $B_z$ data are obtained after portioning the hyperpolarized solution into two \SI{1.0}{\centi\meter\cubed} aliquots, each stored at 5 mT and then placed sequentially next to the zero-field OPM.  Red curves show monoexponential decays with fitted decay time constant $T_1(\mathrm{zero\, field}) = \SI{24+-1}{\second}$.  The lower panel data were obtained with a single aliquot, by alternating 1-s acquisition periods at zero field with 8-s storage periods at \SI{30}{\micro\tesla}. Signal loss in between acquisition periods is attributed to natural decay, fitted by the red curve to $T_1(\SI{30}{\micro\tesla}) = \SI{24+-1}{\second}$, which within error is equal to $T_1(\mathrm{zero\, field})$.
\label{fig:figure3}
\end{figure*}

To apply XY4 to a system of hyperpolarized \textsuperscript{13}C spins, the duration of each dc pulse is increased by a factor $\gamma_{\rm H}/\gamma_{\rm C}\approx 3.97$ to produce rotation angles of $\pi$ and $\sim$$4\pi$ on \textsuperscript{13}C and \textsuperscript{1}H, respectively.  In addition to dynamically decoupling small residual magnetic fields, this combination of angles also decouples the \textsuperscript{13}C-\textsuperscript{1}H scalar ($J_{\rm{CH}}$) coupling interaction when performed sufficiently quickly; $\tau\ll |1/J_{\rm CH}|$. The polarimetry technique can therefore be applied across a wide range of organic molecules.

To demonstrate application to substances of clinical and pre-clinical relevance, we measured the OPM signal during \textsuperscript{13}C XY4 dynamical decoupling for solutions of hyperpolarized [1-\textsuperscript{13}C]-pyruvic acid (PA = \textsuperscript{12}CH\textsubscript{3}{\textsuperscript{12}C}O\textsuperscript{13}COOH; $^3J_{\rm CH}=1.3$ Hz). The samples were prepared initially as \SI{24}{\micro\liter} volumes of 0.2 wt.$\%$ OX-063 radical (GE Healthcare) and 0.09 wt.$\%$ Dotarem (Guerbet S.A., France) dissolved in neat-liquid PA (Sigma Aldrich).  Each of these were polarized at \SI{1.4}{\kelvin} in a magnetic field of \SI{3.5}{\tesla}, using a commercial dissolution-DNP polarizer (HyperSense, Oxford Instruments Ltd.) via 40 minutes of positive-lobe microwave irradiation near 95 GHz.  After polarization, the cryogenic samples were dissolved in a superheated phosphate-buffered-saline solution (4.5 \si{\centi\meter\cubed} in \textsuperscript{1}H\textsubscript{2}O, \SI{10}{\milli\molar} Na\textsubscript{2}HPO\textsubscript{4}, \SI{2}{\milli\molar} KH\textsubscript{2}PO\textsubscript{4}, 1 wt.$\%$ HEPES, 0.01 wt.$\%$ EDTA, 0.1 wt.$\%$ NaCl, 0.2 wt.$\%$ NaOH) and ejected from the instrument at a temperature of \SI{37}{\celsius} into a vial, at a final concentration of $C=0.080(1)$ \si{\molar}.

Immediately after the dissolution process, two \SI{1.0}{\centi\meter\cubed} aliquots of the polarized \textsuperscript{13}C sample were collected and stored in the approximately \SI{5}{mT} fringe magnetic field of the DNP polarizer.  The procedure simulates a situation where fractions of a hyperpolarized sample are inserted into an MRI scanner one after the other\cite{Lees2020NMRBiomedicine}.  The aliquots were then transferred sequentially by syringe into the magnetic shield located approximately 1 m from the polarizer and containing the OPM and sample chamber used for the 2-T-prepolarized water sample.  A solenoid coil provided a guiding magnetic field\cite{Savukov2005PRL} of \SI{30}{\micro\tesla} over the syringe during transport.  For the first aliquot, the guiding field was switched off around 10 seconds after dissolution and the XY4 sequence ($\tau=50$ ms) was initiated and sustained for a time of 6-7 seconds. The liquid was then removed from the apparatus by syringe, the guiding field was switched on again, and the process repeated for the second aliquot; see \autoref{fig:figure3}, upper panel.  During the time at zero field, the OPM signal was recorded, also shown in the upper panel of \autoref{fig:figure3}.  Signals were fit by a monoexponential decay envelope with time constant $T_1 = 24(1)$ s.  The time constant is remarkably long, considering that the paramagnetic Dotarem (Gd\textsuperscript{3+}) and OX-063 species were not removed from the solution. This indicates a low seriousness of paramagnetism-induced relaxation in the zero-to-ultralow-field regime\cite{Bodenstedt2021natcomm}.
One explanation could be that although Gd\textsuperscript{3+} ions remains in the liquid after dissolution, their contact with \textsuperscript{13}C spins is minimized due to chelation in the buffered solution.  The initial amplitude of the \textsuperscript{13}C signal, furthermore, is used to provide a lower-bound value for the polarization.  Amplitudes of the signals for the first and second pyruvate aliquots, 1020 pT and 600 pT respectively, at zero decoupling time equate via \autoref{eq:PA} to \textsuperscript{13}C polarizations of $P=$ 12(1) \% and 7(1) \%, using the 190 pT reference signal amplitude for the 2-T-prepolarized water sample (\autoref{fig:figure2}) and the other parameters as given above. The first one is
in agreement with the \textit{ex post} value $P=$ 13\% calculated from the hyperpolarized vs.\ thermal \textsuperscript{13}C-NMR spectra of a third 1.0 \si{\centi\meter\cubed} aliquot from the same batch, which was also injected at 10.0 s into a bench-top 1.5 T NMR spectrometer (Pulsar, Oxford Instruments Ltd.).

Sample polarimetry is not limited to one-off measurement.  As the lower part of \autoref{fig:figure3} shows, hyperpolarized spins may be sustained across multiple zero-field measurements performed sequentially.  After sample injection into the magnetic shield, the field is repeatedly and rapidly ($<$10 \si{\micro\second}) toggled between the effective-zero field and a \SI{30}{\micro\tesla} holding field to emulate polarization storage in ambient magnetic field.  In the case shown, 8-s storage intervals occur between 1-s zero-field observation windows.  For [\textsuperscript{13}C]-pyruvate, signal decay occurring between windows can be fit to the monoexponential time constant $T_1 = 24$ s, which within error is the same as at zero field under XY4 decoupling.
The observation that the \textsuperscript{13}C $T_1$ is unchanged across three orders of magnitude in field is not so surprising for such a small molecule, given that paramagnetic relaxivities in aqueous solution are generally found to be constant below \SI{1}{\milli\tesla} \cite{Bodenstedt2021natcomm}.  We expect that the use of nonpersistent-radical polarizing agents\cite{Eichhorn2013PNAS110,Pinon2020CommChem3,Gaunt2021ANIE,Negroni2022JACS}, filtration\cite{Gajan2014PNAS} or other strategies to remove paramagnetic dopants from the solution should prolong the decay time at both the storage and detection fields.

Given the growing interest in hyperpolarized MRI, a wide scope of application for the method is envisioned. Solid hyperpolarized [\textsuperscript{13}C]-pyruvate can be transported far away from where it is produced \cite{Ji2017,Capozzi2021CommChem}, for instance to an MRI facility at a hospital located hundreds of kilometers away from the polarizer\cite{Capozzi2022}.
Nondestructive polarization measurement at the use end using a combined dissolution-and-magnetometry device may form part of the quality control procedure carried out upon liquification.
Unlike conventional NMR, our technique does not rely upon maintaining a homogeneous magnetic field around the sample; this means that mandatory procedures for human in-vivo application including sample concentration---via absorption spectroscopy---and pH checks, could be measured in parallel.  This may significantly reduce the delivery time into a human body, which currently stands around 1 minute after dissolution and incurs a substantial natural decay of the polarization\cite{Kurhanewicz2019Neoplasia}.
Hyperpolarization procedures involving chemical reactions of parahydrogen at \si{\micro\tesla} fields\cite{Knecht2021PNAS,Marshall2022arxiv} should also be compatible with detection via OPMs, enabling real-time in-situ monitoring (or optimization) of polarization buildup.

Most liquid-state hyperpolarized NMR and MRI procedures focus on the study of specific chemical sites in isotopically pure compounds (for instance, the terminal carboxyl carbon in pure [1-\textsuperscript{13}C]-PA).  However, it is possible, although uncommon, to hyperpolarize and study more than one isotopomer or compound at the same time, e.g., a solution of [1-\textsuperscript{13}C]-pyruvate and [1-\textsuperscript{13}C]-butyrate\cite{Bastiaansen2016SciRep} or [\textsuperscript{13}C]-urea\cite{Qin2021MRM87}.  The present version of our experiment does not allow polarimetry of two alike spin species at the same time and therefore the magnetometry procedure would have to be performed on the component pure solutions before mixing.  In principle, isotopomer-selective $\pi$ pulses that exploit the chemically specific Zeeman shifts and/or spin-spin J couplings could be used to achieve selectivity in mixtures.  Various other schemes may also be possible, such as magnetic gradiometry of two chemically identical samples that have opposite polarization in one of the compounds.

As a final remark, the square-wave NMR signal may in principle be recorded at any magnetic field, given a sufficiently sensitive detector.  In this work, the nT bias field was chosen to use a zero-field OPM of sensitivity around \SI{10}{\femto\tesla\per\sqrt\hertz} \cite{Tayler2022PRApplied}, although the technical noise background (such as mains hum noise) is much larger, on the order of \SI{1}{\pico\tesla\per\sqrt\hertz}.  Miniaturized versions of zero-field OPMs are produced on a medium commercial scale\cite{quspinQZFMGen3,twinleaf_microserf} for research use in magnetoencephalography\cite{tierney2019optically,Brookes2022TrendsNeurosci} and other biomagnetic imaging studies.  While these must operate in shielded environments for highest performance\cite{Savukov2020JMR,Blanchard2020JMR}, other types of OPMs and magnetic gradiometers can reach sensitivities of tens of \si{\femto\tesla\per\sqrt\hertz} at ambient field\cite{Limes2020} while also being less susceptible to background noise.  In the future, these may remove the current experimental requirement of a magnetic shield.

\begin{figure*}[b]%
\fbox{\includegraphics[]{./toc.png}}%
\caption{Graphical abstract. For Table Of Contents only.
\label{fig:toc}
\end{figure*}

\section*{Acknowledgments}
The work described is funded by:
EU H2020 Marie Sk{\l}odowska-Curie Actions project ITN ZULF-NMR  (Grant Agreement No. 766402); the European Union’s Horizon 2020 research and innovation program (GA-863037);
the Spanish Ministry of Science MCIN with funding from European Union NextGenerationEU (PRTR-C17.I1) and by Generalitat de Catalunya ``Severo Ochoa'' Center of Excellence CEX2019-000910-S; projects OCARINA (PGC2018-097056-B-I00), SAPONARIA (PID2021-123813NB-I00) and MARICHAS (PID2021-126059OA-I00),  funded by MCIN/ AEI /10.13039/501100011033/ FEDER ``A way to make Europe'';
Generalitat de Catalunya through the CERCA program;
Ag\`{e}ncia de Gesti\'{o} d'Ajuts Universitaris i de Recerca Grant Nos. 2017-SGR-1354 and 2021 FI\_B\_01039;
Secretaria d'Universitats i Recerca del Departament d'Empresa i Coneixement de la Generalitat de Catalunya, co-funded by the European Union Regional Development Fund within the ERDF Operational Program of Catalunya (project QuantumCat, ref. 001-P-001644);
Fundaci\'{o} Privada Cellex;
Fundaci\'{o} Mir-Puig;
and the BIST--“la Caixa” initiative in Chemical Biology (CHEMBIO).
MCD Tayler and I Marco-Rius acknowledge financial support through the Junior Leader Postdoctoral Fellowship Programme from ``La Caixa'' Banking Foundation (projects LCF/BQ/PI19/11690021 and LCF/BQ/ PI18/11630020, respectively).

\section*{Supporting Information}
\noindent Additional experimental details including photographs of the experimental setup (PDF).

\section*{Competing interests}
\noindent The authors declare no competing interests.

\section*{References}


\section{Introduction}
\label{section:introduction}

How the cold, empty universe after inflation has evolved to our observable universe has been the subject of rigorous investigation in the past. Aside from reheating \cite{Linde:1990flp}, where the inflaton decays perturbatively to other forms of matter, it is now known that coupling of another field to the inflaton, while it oscillates around its minimum, can lead to a burst of particle production. This is known as parametric resonance or stochastic preheating \cite{Traschen:1990sw, Shtanov:1994ce, Kofman:1994rk, Kofman:1997yn}. During this period, the number density of the preheat fields can grow or decay due to the expansion of the universe, although it mainly grows most of the time. The products of this process finally interact with each other and thermalize at a temperature known as reheating temperature. In most inflationary models, even the ones that have roots in string theory, the existence and coupling of such preheat fields are assumed in an ad-hoc manner.

An inflationary model which has well-based theoretical motivations from string theory is Matrix inflation, or simply M-flation \cite{Ashoorioon:2009wa, Ashoorioon:2011ki}. Different aspects of this model were investigated in the literature \cite{Ashoorioon:2009sr, Ashoorioon:2013oha, Ashoorioon:2014jja}. In this class of models, we deal with three $N\times N$ Hermitian matrices $\Phi_{i}$ ($i=1,2,3$), and the potential is assumed to be a function of $\Phi_{i}$ or the commutations $\left[\Phi_{i},\Phi_{j}\right]$. Therefore, the potential of $\Phi_{i}$ appears in these models as a function of $\mathrm{Tr}\,\Phi_{i}^{2}$, $\mathrm{Tr}\,\epsilon_{ijk}\Phi_{i}\left[\Phi_{j},\Phi_{k}\right]$, and $\mathrm{Tr}\,\left[\Phi_{i},\Phi_{j}\right]^{2}$. In Ref. \cite{Ashoorioon:2009wa, Ashoorioon:2011ki} it was discussed how the potential can be motivated from the brane dynamics in string theory.

The M-flation scenario can resolve the fine-tuning problem of the self-couplings in chaotic inflation models \cite{Linde:1983gd}, as it has been shown explicitly in \cite{Ashoorioon:2009wa}. In the setting of standard chaotic inflation, in order to make the simple chaotic inflation models compatible with the CMB observations, the potential parameters are required to be very small which is not acceptable from the point of view of quantum field theory. In addition, the inflaton field displacement in the chaotic setups is super-Planckian which is not favored by the quantum field theory considerations, and it also is not consistent with the swampland distance conjecture \cite{Ooguri:2018wrx, Agrawal:2018own}. The M-flation scenario, however, can resolve these problems by providing the potential parameters of order unity, and also presenting sub-Planckian field displacements during inflation \cite{Ashoorioon:2009wa}. In this setup, like other multi-field inflation models, the isocurvature perturbations are present along with the adiabatic one. Of course, due to their masses, their amplitudes decay exponentially by the end of inflation. Since the coupling of these isocurvature fields is related to the self-couplings of the inflaton, their amplitude can be understood exactly. Also at the end of inflation, when the slow-roll condition is violated, these isocurvature fields can act as preheat fields and drain energy from the inflation. Knowing their coupling to the inflaton, one can have exact predictions for the gravitational wave spectrum produced during preheating \cite{Ashoorioon:2013oha}.

Despite these achievements, M-flation suffers from several shortcomings. For instance, although this model is based on string theory, so far the issue of compactification to four dimensions has not been worked out. More strictly, thus far, it is not clear how the 4-dimensional gravity arises when one compactifies the extra dimensions. In addition, in order to solve the problems of the smallness of the potential self-couplings from the bare ones, the ranks of the matrices, which represent the number of D3-branes in the stack, are required to be of order $N\sim10^{4}-10^{5}$, which is large and makes the backreaction on the underlying background geometry an important issue that needs to be addressed. It is also not compatible with the latest bound on the tensor-to-scalar ratio by the Planck data \cite{Planck:2018jri}.

In the setup of non-minimal M-flation (non-$\mathbb{M}$-flation) \cite{Ashoorioon:2019kcy}, in which the inflaton field couples non-minimally to gravity, we can resolve some of these problems. The non-minimal coupling to gravity is something that has been observed in the renormalization of the gravity in presence of matter. In fact, it has been shown that the loop correction to the graviton-scalar-scalar vertex generates a term proportional to $\ensuremath{\xi\Lambda^{2}\phi^{2}/M_P^{2}}$ with $\xi\simeq \mathcal{O}(1)$ \cite{Ashoorioon:2011aa}, where $\Lambda$ is the cutoff of the theory. Noting that $\Lambda$ is at most $M_P$, the resulting coupling is of order one at most, and cannot justify the large values of $\xi$ if it is needed. In $\mathbb{K}$L$\mathbb{M}$T \cite{Kachru:2003sx} though, when the superpotential is proportional to the volume modulus, it was possible to obtain large non-minimal coupling parameters \cite{Ashoorioon:2019kcy}. One should note that the idea of the non-minimal coupling to the gravitational field has been adopted in the Higgs inflation model \cite{Bezrukov:2007ep}. It will be more convenient to analyze the model, if we go through a conformal transformation to the Einstein frame. In non-$\mathbb{M}$-flation, it is possible to improve the observational consistency of the symmetry-breaking potential with the Planck 2018 CMB data \cite{Planck:2018jri} compared to the minimal M-flation predictions \cite{Ashoorioon:2019kcy}. The number of the D3-branes in non-$\mathbb{M}$-flation is obtained to be of order $N\lesssim10^{2}$ which is considerably less than the required number in the minimal M-flation setup, which is one appealing consequence of non-minimal coupling in the model \cite{Ashoorioon:2019kcy}. Another appealing feature of the non-$\mathbb{M}$-flation, in the supersymmetric setup, was the hope that preheating around the symmetry-breaking vacuum could be addressed. However as we will revisit this problem here, we notice that the problem persists in the non-minimal symmetry-breaking M-flation. That is why we turn our attention to a supersymmetric setup in this paper.

As noted above, after inflation, the spectator fields may play an important role in the Universe's evolution. In particular, their coupling to the inflaton field can provide a successful mechanism of parametric resonance that explains the particle production in the post-inflationary universe. This process is also known as the preheating process during which the universe experiences a decelerating expansion \cite{Kofman:1994rk, Kofman:1997yn}. During preheating, the number of particles can increase significantly. The particles produced at this stage are far away from thermal equilibrium and have enormously large occupation numbers. In the next stage, the interactions of the existing fields with the other fields become important, and the previously produced particles decay. In the final step, the cosmic fluid will be thermalized, and this process generally takes a very long time which is inversely proportional to the coupling constants.

The quantum dispersion of the preheat field squared leads to the modification of the potential which changes the evolution of the post-inflationary universe. The phenomenon was mainly noticed in the symmetry-breaking potentials where the quantum fluctuations of the preheat field restore the symmetry and lead to the production of topological defects and secondary stage of inflation \cite{Kofman:1995fi, Felder:2000sf}. The period of accelerated expansion which takes place after preheating may have significant theoretical and observational implications. In particular, such a transient period of acceleration could be crucial in solving the moduli and gravitino problems. This phase of acceleration can be regarded as secondary non-thermal inflation and it may be an alternative to the thermal inflation \cite{Lyth:1995hj} proposed for solving these problems \cite{Felder:2000sf}. In the example that we consider here, a similar phenomenon is noted for a light massive preheat field that acquires a dispersion during inflation preheating. The mass term acts for a period of time like a cosmological constant that lifts the inflaton potential and leads to accelerated expansion.

Preheating after inflation increases the number density of particles in some specific momentum instability bands. This corresponds to inhomogeneities that source gravitational waves at the second order in perturbation theory \cite{Easther:2006vd, Ashoorioon:2022raz}. The spectrum of the stochastic background of the gravitational wave is different from those produced by a first-order phase transition \cite{Ashoorioon:2022raz}. In \cite{Easther:2006vd}, the GWs from preheating have been investigated in a simple inflationary setup with the quadratic potential $V(\phi)=\mu^{2}\phi^{2}/2$ where $\mu$ is a constant parameter. In their investigation, they have shown that for different values of the inflaton effective mass in the wide range $\mu = 10^{-18}-10^{-6}m_{P}$ where $m_P$ is the Planck mass, the amplitude of the GWs is almost identical and it is of order $\Omega_{\mathrm{GW}}h^{2}\sim10^{-11}-10^{-10}$, but the spectrum appears in different frequencies. For the important case of $\mu=10^{-6}m_{P}$ the spectrum appears in the frequencies $f\sim10^{8}-10^{10}\mathrm{Hz}$. In this paper, we refer to this case as the benchmark reheating scenario, and we will compare our findings with the results of that model\footnote{In our recent paper \cite{Ashoorioon:2022raz}, it was shown that in the lower energy scales, the amplitude of the GWs spectrum shows some dependency on the parameter $\mu$.}.

We study the effect of preheating process on the inflaton potential in the framework of non-$\mathbb{M}$-flation, but away from the symmetry-breaking configuration. The potential thus has only one minimum. We use the numerical code LATTICEEASY \cite{Felder:2000hq} to compute the root mean square (r.m.s.) of the scalar spectator field in our model, and then we use it to calculate the modified term in the potential. We include the contribution of the additional potential energy in the ensuing background dynamics. Depending on the potential parameters and the properties of the dominant preheat field, there will be different preheating scenarios in our setup. We show that in a specific case, the lifted potential is capable of providing a transition to a temporary accelerating phase of expansion. In this paper, we show that this transient period of accelerated expansion, suppresses the number density of the produced preheat particles and hence reduces the amplitude of the gravitational waves generated during preheating.

This paper is structured as follows. First, in Sec. \ref{section:model}, we review briefly the basics of the non-$\mathbb{M}$-flation model. Next, in Sec. \ref{section:preheating}, we present the basic equations for the study of preheating in this setup. Subsequently in Sec. \ref{section:transient_acceleration}, we will consider a special case in the setting of non-$\mathbb{M}$-flation that is capable of providing a transient acceleration after preheating. We will then investigate the effect of the transient acceleration on the number density and also the spectrum of the stochastic background of GWs generated during preheating. Finally, we summarize our findings in Sec. \ref{section:conclusions}.

\section{Non-$\mathbb{M}$-flation}
\label{section:model}

In non-$\mathbb{M}$-flation, it is assumed that the matrix fields $\Phi_{i}$ couple non-minimally to the gravitational field, so the action of the scalar fields of the model in the Jordan frame can be written as
\begin{equation}
\label{SJ-Phi}
S_{J}=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}\left(1+\xi\sum\limits _{i=1}^{3}\mathrm{Tr}\left(\Phi_{i}^{2}\right)\right)R-\frac{1}{2}\sum\limits _{i=1}^{3}{\rm Tr}\left(\partial_{\mu}\Phi_{i}\partial^{\mu}\Phi_{i}\right)-V\left(\Phi_{i},\left[\Phi_{i},\Phi_{j}\right]\right)\right] \; .
\end{equation}
Throughout this paper, we work in the units where the reduced Planck mass is set to unity, $M_{P}\equiv1/\sqrt{8\pi G}=1$. The coupling strength of the matrix fields with gravitation is determined by the coupling constant $\xi$. We do not include the gauge fields in the present analysis. The inflationary potential $V$ in this action is a function of the matrices $\Phi_{i}$, and their commutations $\left[\Phi_{i},\Phi_{j}\right]$. In Ref. \cite{Ashoorioon:2009wa} it has been discussed that it is possible to derive the potential from the examination of the string theory branes. In that reference, it has also shown that in the special context of string theory and to the dominant order of $\Phi_{i}$ and $\left[\Phi_{i},\Phi_{j}\right]$, the potential takes the following form
\begin{equation}
\label{V-Phii}
V\left(\Phi_{i},\left[\Phi_{i},\Phi_{j}\right]\right)=\mathrm{Tr}\left(-\frac{\lambda}{4}\left[\Phi_{i},\Phi_{j}\right]\left[\Phi_{i},\Phi_{j}\right]+\frac{i\kappa}{3}\epsilon_{jkl}\left[\Phi_{k},\Phi_{l}\right]\Phi_{j}+\frac{m^{2}}{2}\Phi_{i}^{2}\right) \; ,
\end{equation}
where $\lambda$ is a dimensionless constant, while the constants $\kappa$ and $m$ have dimensions of mass. These constants have some interpretations from string theory: the constant $\lambda=8\pi g_{s}=2g_{YM}^{2}$ is related to the string coupling constant $g_s$ or the Yang-Mills constant $g_{YM}$, the constant $\kappa=\hat{\kappa}g_{s}\sqrt{8\pi g_{s}}$ has relation with the Ramond-Ramond antisymmetric form strength $\hat{\kappa}$, and the constant $m$ comes from the three spatial coordinates along the D3-branes in the metric of the background SUGRA theory \cite{Ashoorioon:2009wa, Ashoorioon:2011ki}. The background is a solution to SUGRA equations with a constant dilaton if $\lambda m^{2}=4\kappa^{2}/9$. However in this paper, we consider more general configurations and deviate from this condition.

Since the matrices $\Phi_{i}$s are $N-$dimensional, we deal with $3N^{2}$ real scalar fields, and hence the study of the model in its most general form becomes very complicated. Instead, we introduce the matrices $J_{i}$ ($i=1,2,3$) that are the $N\times N$ generators of the SU(2) algebra, and satisfy the commutation relation $\left[J_{i},J_{j}\right]=i\,\epsilon_{ijk}\,J_{k}$. Thus, we can decompose the inflaton matrices into two parts,
\begin{equation}
\label{Phii-phihat}
\Phi_{i}=\hat{\phi}J_{i}+\Psi_{i} \, .
\end{equation}
The first and second parts in this decomposition are respectively the parallel and perpendicular components to the $N\times N$ matrix representation of the SU(2) algebra (i.e. $\mathrm{Tr}\left(J_{i}\Psi_{i}\right)=0$). In Refs. \cite{Ashoorioon:2009wa, Ashoorioon:2013oha}, it has been argued that if the fields $\Psi_{i}$ are turned off at the early stages, they remain classically zero during the inflation. Thus we call them the spectator fields. Ignoring the spectator fields, the inflationary trajectory is determined by the $\hat{\phi}$ field which is the component of the matrices along the SU(2) generator. Accordingly, the action \eqref{SJ-Phi} can be rewritten as
\begin{equation}
\label{SJ-phihat}
S_{J}=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}\left(1+\frac{\xi}{M_{P}^{2}}\mathrm{Tr}\left(J_{i}^{2}\right)\hat{\phi}^{2}\right)R+\mathrm{Tr}\left(J_{i}^{2}\right)\left(\frac{1}{2}\left(\frac{d\hat{\phi}}{dt}\right)^{2}-\frac{\lambda}{2}\hat{\phi}^{4}+\frac{2\kappa}{3}\hat{\phi}^{3}-\frac{m^{2}}{2}\hat{\phi}^{2}\right)\right] \, ,
\end{equation}
with ${\rm Tr}\left(J_{i}^{2}\right)=N\left(N^{2}-1\right)/4$. Since the kinetic term in this action is non-canonical, we redefine the scalar field as
\begin{equation}
\label{phi-phihat}
\phi=\sqrt{\mathrm{Tr}\left(J_{i}^{2}\right)}\,\hat{\phi} \, .
\end{equation}
As a result, the action takes the following form
\begin{equation}
\label{SJ}
S_{J}=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}\left(1+\xi\phi^{2}\right)R+\frac{1}{2}\left(\frac{d\phi}{dt}\right)^{2}-V_{0}(\phi)\right] \, .
\end{equation}
Introducing the effective coupling constants $\lambda_{\mathrm{eff}}\equiv2\lambda/\mathrm{Tr}\left(J_{i}^{2}\right)$ and $\kappa_{\mathrm{eff}}\equiv\kappa/\sqrt{\mathrm{Tr}\left(J_{i}^{2}\right)}$, the inflaton potential \eqref{V-Phii} can be written as
\begin{equation}
\label{V0}
V_{0}(\phi)=\frac{\lambda_{\mathrm{eff}}}{4}\phi^{4}-\frac{2\kappa_{\mathrm{eff}}}{3}\phi^{3}+\frac{m^{2}}{2}\phi^{2} \, .
\end{equation}
To get rid of the non-minimal coupling to gravity, following Ref. \cite{Bezrukov:2007ep}, we do the following conformal transformation from the Jordan to the Einstein frame,
\begin{equation}
\label{gtildemunu}
\tilde{g}_{\mu\nu}=\Omega^{2}g_{\mu\nu}\, , \qquad \Omega^{2}=1+\xi\phi^{2} \, .
\end{equation}
The tilde represents the quantities in the Einstein frame. This conformal transformation makes the kinetic term in the Jordan frame to be non-canonical, and to resolve this problem, we introduce the new scalar field $\chi$ which is the scalar field in the Einstein frame, and it has a relation to the Jordan frame scalar field $\phi$ as
\begin{equation}
\label{dchidphi}
\frac{d\chi}{d\phi}=\frac{\sqrt{\Omega^{2}+6\xi^{2}\phi^{2}}}{\Omega^{2}}=\frac{\sqrt{\xi\phi^{2}(6\xi+1)+1}}{\xi\phi^{2}+1} \, .
\end{equation}
Therefore, the action in the Einstein frame turns into
\begin{equation}
\label{SE}
S_{E}=\int d^{4}x\sqrt{-\tilde{g}}\left[\frac{1}{2}\tilde{R}+\frac{1}{2}\left(\frac{d\chi}{dt}\right)^{2}-U(\chi)\right] \, ,
\end{equation}
where $U(\chi)$ is the inflationary potential in the Einstein frame and for it we have
\begin{equation}
\label{U}
U(\chi)=\frac{V_{0}\left(\phi(\chi)\right)}{\Omega^{4}\left(\phi(\chi)\right)} \, .
\end{equation}
From the ordinary differential equation \eqref{dchidphi}, the field $\chi$ can be expressed in terms of the field $\phi$ as
\begin{equation}
\label{chi-phi}
\chi=f(\phi)=\sqrt{\frac{1}{\xi}+6}\,\sinh^{-1}\left[\sqrt{\xi\left(6\xi+1\right)}\,\phi\right]-\sqrt{6}\,\tanh^{-1}\left[\frac{\sqrt{6}\,\xi\text{\ensuremath{\phi}}}{\sqrt{\xi(6\xi+1)\text{\ensuremath{\phi}}^{2}+1}}\right] \, .
\end{equation}
It is not simple to solve the above equation to express $\phi$ in terms of $\chi$, and for this purpose, we can use the inverse function,
\begin{equation}
\label{phi-chi}
\phi=f^{-1}(\chi) \, .
\end{equation}
In the numerical calculations, it is more convenient to use the following function which provides an appropriate approximation for $\phi$ versus $\chi$
\begin{equation}
\label{phi-chi-approx}
\phi\approx\frac{a_{1}\chi+a_{2}\chi^{2}+a_{3}\chi^{3}+a_{4}\chi^{4}+a_{5}\chi^{5}+a_{6}\chi^{6}+a_{7}\chi^{7}+a_{8}\chi^{8}+a_{9}\chi^{9}+a_{10}\chi^{10}}{1+b_{1}\chi+b_{2}\chi^{2}+b_{3}\chi^{3}+b_{4}\chi^{4}+b_{5}\chi^{5}+b_{6}\chi^{6}+b_{7}\chi^{7}+b_{8}\chi^{8}+b_{9}\chi^{9}+b_{10}\chi^{10}} \, .
\end{equation}
The numerical coefficients in this equation can be determined by curve-fitting, for each set of model parameters.

\section{Preheating in non-$\mathbb{M}$-flation}
\label{section:preheating}

In this section, we aim to study preheating in the framework of non-$\mathbb{M}$-flation. For this purpose, we take the spectator scalar fields into account, and then the Jordan frame action will be
\begin{equation}
\label{SJ-Psi}
S_{J}=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}\left(1+\xi\phi^{2}\right)R-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V_{0}(\phi)-\frac{1}{2}\sum\limits _{i}g^{\mu\nu}\partial_{\mu}\Psi_{i}\partial_{\nu}\Psi_{i}-V_{(2)}\left(\phi,\Psi_{i}\right)\right] \, .
\end{equation}
The potential in this equation has been decomposed into two parts,
\begin{equation}
\label{V}
V=V_{0}\left(\phi\right)+V_{(2)}\left(\phi,\Psi_{i}\right) \, .
\end{equation}
The potential $V_{0}\left(\phi\right)$ has already given by Eq. \eqref{V0}, and for the potential $V_{(2)}\left(\phi,\Psi_{i}\right)$, we have
\begin{equation}
\label{V2}
V_{(2)}\left(\phi,\Psi_{i}\right)=\frac{1}{2}M_{\Psi}^{2}(\phi)\Psi_{i}^{2} \, ,
\end{equation}
where
\begin{equation}
\label{MPsi2}
M_{\Psi}^{2}(\phi)=\frac{1}{2}\lambda_{{\rm eff}}\omega(\omega-1)\phi^{2}+2\kappa_{{\rm eff}}\omega\phi+m^{2} \, .
\end{equation}
The parameter $\omega$ takes integers values, and we can calculate $M_\Psi$ which is the effective mass of the $\Psi_i$ modes. There are two sets of solutions for the effective mass of the scalar spectator fields that are known as $\alpha$-modes and $\beta$-modes. For the $\alpha$-modes we have $\omega=-(j+2)$ where $j$ is an integer in the range $0\leq j\leq N-2$. The degeneracy of each $\alpha_j$-mode is $2j+1$, and so the total number of these modes is equal to $(N-1)^2$ \footnote{The are $N^2-1$ zero modes which are replaced with the gauge modes too. We do not take them into account in this work.}. The $\alpha$-mode with $j = 0$ corresponds to the quantum fluctuations in the SU(2) direction and it is an adiabatic mode. By discarding the contribution of this mode, the total number of degeneracy becomes $(N-1)^2 - 1$. For the $\beta$-modes, we have $\omega = j-1$ with $1\leq j\leq N$. Each $\beta$-mode has $2j+1$ degeneracy, and hence the total degeneracy will be resulted in to be $(N+1)^2 - 1$.

Exploiting the conformal transformation \eqref{gtildemunu}, we can simply show
\begin{equation}
\partial_{\mu}\Psi_{i}\partial^{\mu}\Psi_{i}=g^{\mu\nu}\partial_{\mu}\Psi_{i}\partial_{\nu}\Psi_{i}=\Omega^{2}\tilde{g}^{\mu\nu}\partial_{\mu}\Psi_{i}\partial_{\nu}\Psi_{i} \, .
\end{equation}
Substituting this into Eq. \eqref{SJ-Psi}, we obtain the action in the Einstein frame as
\begin{equation}
\label{SE-Psi}
S_{E}=\int d^{4}x\sqrt{-\tilde{g}}\left[\frac{1}{2}\tilde{R}-\frac{1}{2}\tilde{g}^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi-U(\chi)-\frac{1}{2\Omega^{4}}\sum\limits _{i}g^{\mu\nu}\partial_{\mu}\Psi_{i}\partial_{\nu}\Psi_{i}-\tilde{V}_{(2)}\left(\chi,\Psi_{i}\right)\right] \, .
\end{equation}
The new potential in this action will be
\begin{equation}
\label{V2tilde}
\tilde{V}_{(2)}\left(\chi,\Psi_{i}\right)=\frac{V_{(2)}\left(\chi,\Psi_{i}\right)}{\Omega^{4}\left(\phi(\chi)\right)}=\frac{1}{2}\frac{M_{\Psi_{i}}^{2}\left(\phi(\chi)\right)}{\Omega^{4}\left(\phi(\chi)\right)}\Psi_{i}^{2} \, .
\end{equation}
Imposing the redefinition $\tilde{\Psi}_{i}\equiv\Psi_{i}/\Omega$ and also using the relation $\tilde{g}^{\mu\nu}=g^{\mu\nu}/\Omega^{2}$, action \eqref{SE-Psi} turns into
\begin{equation}
\label{SE-Psitilde}
S_{E}=\int d^{4}x\sqrt{-\tilde{g}}\left[\frac{1}{2}\tilde{R}-\frac{1}{2}\tilde{g}^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi-U(\chi)-\frac{1}{2}\sum\limits _{i}\tilde{g}^{\mu\nu}\partial_{\mu}\tilde{\Psi}_{i}\partial_{\nu}\tilde{\Psi}_{i}-\bar{V}_{(2)}\left(\chi,\tilde{\Psi}_{i}\right)+\frac{\dot{\Omega}}{\Omega}\sum\limits _{i}\tilde{\Psi}_{i}\dot{\tilde{\Psi}}_{i}\right] \, .
\end{equation}
The redefined potential in this equation is given by
\begin{equation}
\label{V2bar}
\bar{V}_{(2)}\left(\chi,\tilde{\Psi}_{i}\right)=\frac{1}{2}M_{\tilde{\Psi}}^{2}\tilde{\Psi}_{i}^{2} \, ,
\end{equation}
where
\begin{equation}
\label{MPsit2}
M_{\tilde{\Psi}}^{2}=\frac{1}{\Omega^{2}}\left(M_{\Psi}^{2}-\dot{\Omega}^{2}\right) \, .
\end{equation}
In the above equation, if the second term in the parentheses is much less than the first term, $\dot{\Omega}^{2}\ll M_{\Psi}^{2}$, we will have
\begin{equation}
\label{MPsit2-approx}
M_{\tilde{\Psi}}^{2}\approx\frac{M_{\Psi}^{2}}{\Omega^{2}} \, .
\end{equation}
The use of this approximation will simplify our calculations considerably, in particular for the lattice simulations. We will use this approximation in the case we analyze, and then in the subsequent steps, we verify that it is indeed valid.

The mode $\tilde{\Psi}_{i}(t,\mathbf{x})$ can be decomposed into its Fourier components as follows
\begin{equation}
\label{Psiti-Psitik}
\tilde{\Psi}_{i}(t,\mathbf{x})=\int\frac{d^{3}k}{(2\pi)^{3/2}}\left[\tilde{\Psi}_{i_{k}}(t)\hat{a}_{k}e^{-i\mathbf{k.x}}+\tilde{\Psi}_{i_{k}}^{*}(t)\hat{a}_{k}^{\dagger}e^{i\mathbf{k.x}}\right] \, ,
\end{equation}
where $\hat{a}_{k}$ and $\hat{a}_{k}^{\dagger}$ are the annihilation and creation quantum operators, respectively. In addition, we rescale the Fourier mode of the spectator field as $\bar{\Psi}_{i_{k}}\equiv a^{3/2}\tilde{\Psi}_{i_{k}}$, and consequently arrive at the following evolution equation \cite{Ashoorioon:2019kcy}
\begin{equation}
\label{ddotPsibik}
\ddot{\bar{\Psi}}_{i_{k}}+\omega_{k}^{2}\bar{\Psi}_{i_{k}}=0 \, .
\end{equation}
The time-dependent angular frequency in this equation is given by
\begin{equation}
\label{omegak}
\omega_{k}^{2}\equiv\frac{k^{2}}{a^{2}}+\frac{M_{\Psi}^{2}}{\Omega^{2}}-\frac{3}{4}H^{2}-\frac{3}{2}\frac{\ddot{a}}{a}+3H\frac{\dot{\Omega}}{\Omega}-2\left(\frac{\dot{\Omega}}{\Omega}\right)^{2}+\frac{\ddot{\Omega}}{\Omega} \, .
\end{equation}
To solve the second-order different equation \eqref{ddotPsibik}, we impose the Bunch-Davies vacuum initial condition on the rescaled spectator mode at the start of preheating,
\begin{equation}
\label{PsibikBD}
\bar{\Psi}_{i_{k}}\to\frac{1}{\sqrt{2\omega_{k}}}e^{-i\omega_{k}t} \, .
\end{equation}
The solution of the differential equation \eqref{ddotPsibik} is usually used to calculate the number density of the produced particles which is given by
\begin{equation}
\label{nk}
n_{k}=\frac{\omega_{k}}{2}\left(\frac{1}{\omega_{k}^{2}}\left|\dot{\bar{\Psi}}_{i_{k}}\right|^{2}+\left|\bar{\Psi}_{i_{k}}\right|^{2}\right)-\frac{1}{2} \, .
\end{equation}
To solve the differential equation \eqref{ddotPsibik}, the value of the comoving wavenumber $k$ should be determined, and in our work, we take it as $k = 0$. The number density is an important criterion that implies the extent of the efficiency of the preheating process.

Due to quantum fluctuations of the scalar fields during preheating the effective potential acquires the following correction \cite{Kofman:1995fi, Felder:2000sf}
\begin{equation}
\label{DeltaU}
\Delta U=\frac{1}{2}M_{\tilde{\Psi}}^{2}\left\langle \tilde{\Psi}_{i}^{2}\right\rangle  \, .
\end{equation}
In our work, to calculate the dispersion of the spectator field, $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle$, we perform a lattice simulation by using the computational code LATTICEEASY \cite{Felder:2000hq}. Including the above correction term, the effective potential will take the following form
\begin{equation}
\label{Utot}
U_{\mathrm{tot}}(\chi)=U(\chi)+\Delta U \, .
\end{equation}
The correction term in the potential energy may have a considerable impact on the dynamics of the universe during preheating. This term may even cause the universe undergoes a transient accelerating expansion during preheating, and this in turn may have cosmological outcomes. The effect of this correction on preheating will be examined in detail in the next section for a special scenario of the non-$\mathbb{M}$-flation model.

\section{Transient acceleration during preheating in non-$\mathbb{M}$-flation}
\label{section:transient_acceleration}

In the previous section, we have presented the basic equations required for the investigation of preheating in the framework of non-$\mathbb{M}$-flation. In this section, we consider a specific case of this setting and use those equations to investigate the preheating process for it.

It should be noted that taking non-vanishing values for the parameter $\kappa_{\mathrm{eff}}$ leads to some problems. If we take this parameter as $\lambda m^{2}=4\kappa^{2}/9$ which restricts the inflaton dynamics to the super-gravity equations, then the model fails to provide a successful process of preheating in the case that the inflaton field never reaches the super-symmetry minimum. This case will be discussed in detail in Appendix \ref{appendix:symmetry-breaking}. In addition, in the cases for which the inflaton can cross the potential bump and reach the super-symmetry minimum, then the effective mass squared of the preheat field will be negative in a range of $\chi$. This in turn cause that after inclusion of the potential correction term in the background dynamics, the effective potential of inflaton acquires negative values in that range of $\chi$. The negative values of the effective potential lead to some critical problems in numerical computations, and therefore here we focus on the case $\kappa_{\mathrm{eff}}=0$ to prevent the numerical errors in our computations.

Also, it should be noted that if we take $\mu\gtrsim0.001M_{P}$, then the dispersion of the Einstein-frame preheat field, $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle $, gets very large values which are even of order $M_P$ or more. This means that the magnitude of quantum fluctuations of the preheat field becomes of the order of the inflaton field magnitude, and this in turn makes some problems in the LATTICEEASY code. Therefore, in order to do a successful lattice simulation, we have taken $\mu=10^{-4}M_{P}$ which is small enough to avoid such problems. For the coupling constant to gravity, we choose $\xi=100$ which is large enough to make the inflationary observables compatible with the current CMB observations. We further take $\omega=-20$ which is related to $\alpha$-mode with $j = 18$.

Since we have set $\kappa_{\mathrm{eff}}=0$ here, then the potential in both the Jordan and Einstein frames has only one minimum. In the Jordan frame, the potential  is obtained from Eq. \eqref{V0} as
\begin{equation}
\label{V0-2}
V_{0}(\phi)=\frac{\text{\ensuremath{\lambda}}_{\mathrm{eff}}}{4}\phi^{2}\left(\phi^{2}+\mu^{2}\right) \, .
\end{equation}
Substituting this into Eq. \eqref{U}, the Einstein-frame potential is obtained as
\begin{equation}
\label{U-2}
U(\chi)=\frac{\text{\ensuremath{\text{\ensuremath{\lambda}}_{\mathrm{eff}}}}\phi^{2}(\chi)\left[\phi^{2}(\chi)+\mu^{2}\right]}{4\left[\xi\phi^{2}(\chi)+1\right]^{2}} \, .
\end{equation}
The potential of this case is quite symmetric around its minimum in both the Jordan and Einstein frames. Because of the symmetry of this potential, both the left and right branches of the potential give rise to similar inflationary scenarios. However,  here we concentrate on the inflationary trajectory from the right branch of the potential. From fixing the amplitude of the scalar perturbations according to the recent CMB observations, we find $\text{\ensuremath{\lambda}}_{\mathrm{eff}}=3.780\times10^{-6}$. This value of $\text{\ensuremath{\lambda}}_{\mathrm{eff}}$ gives the number of simultaneous branes as $N\approx130$ which is much less than what would be required in the M-flation scenario. We further get the scalar spectral index and tensor-to-scalar ratio, respectively, as $n_{s}=0.9678$ and $r=0.0030$ which are in good agreement with the Planck 2018 observations \cite{Planck:2018jri}.

To do the numerical commutations for this case, it is convenient to use the approximated function in Eq. \eqref{phi-chi-approx} to express $\phi$ in terms of $\chi$. The constant coefficients of this function for this case are presented in Table \ref{table:a-b}. Using these coefficients, the plot of $\phi$ against $\chi$ is obtained as depicted in Fig. \ref{figure:phi-chi}. In the figure, we have also presented the result coming from the direct calculation of $\phi=f^{-1}(\chi)$ provided by a numerical root-finding method. The figure implies that function \eqref{phi-chi-approx} provides a very good approximation for $\phi$, and so we can use this function for numerical purposes in the LATTICEEASY code. Also, it is much appropriate to apply approximation \eqref{MPsit2-approx} for $M_{\tilde{\Psi}}^{2}$ in the code. In the following, we will verify explicitly that this approximation is really valid in our examination. In our lattice simulation, we have set the number of points along each edge of the cubical lattice and the size of the box (i.e. length of each edge) in the lattice units as $N = 128$ and $L = 40$, respectively.

\begin{table}[t]
\centering
\scalebox{1.0}{
\begin{tabular}{|c|c|c|c|}
\hline
$a_{1}$ & $9.263\times10^{-1}$ & $b_{1}$ & $-3.343\times10^{-7}$\tabularnewline
\hline
$a_{2}$ & $-2.562\times10^{-7}$ & $b_{2}$ & $8.087\times10^{3}$\tabularnewline
\hline
$a_{3}$ & $4.042\times10^{3}$ & $b_{3}$ & $-5.910\times10^{-4}$\tabularnewline
\hline
$a_{4}$ & $-1.802\times10^{-4}$ & $b_{4}$ & $3.539\times10^{6}$\tabularnewline
\hline
$a_{5}$ & $9.761\times10^{5}$ & $b_{5}$ & $1.367\times10^{-2}$\tabularnewline
\hline
$a_{6}$ & $9.271\times10^{-3}$ & $b_{6}$ & $1.286\times10^{8}$\tabularnewline
\hline
$a_{7}$ & $2.079\times10^{7}$ & $b_{7}$ & $2.072$\tabularnewline
\hline
$a_{8}$ & $3.344\times10^{-1}$ & $b_{8}$ & $4.225\times10^{8}$\tabularnewline
\hline
$a_{9}$ & $4.102\times10^{7}$ & $b_{9}$ & $5.963$\tabularnewline
\hline
$a_{10}$ & $5.548\times10^{-1}$ & $b_{10}$ & $8.054\times10^{6}$\tabularnewline
\hline
\end{tabular}
}
\caption{The constant coefficients in function \eqref{phi-chi-approx} which is used in our investigation to approximate the Jordan-frame scalar field $\phi$ in terms of the one in the Einstein frame $\chi$.}
\label{table:a-b}
\end{table}

\begin{figure}
\centering
\includegraphics[scale=0.4]{phi-chi-2.pdf}
\caption{The scalar field in the Jordan frame ($\phi$) in terms of the one in the Einstein frame ($\chi$) in our setting. The solid curve represents the result from the numerical root-finding method to evaluate $\phi=f^{-1}(\chi)$, but the dashed curve is generated by using the approximate function \eqref{phi-chi-approx}.}
\label{figure:phi-chi}
\end{figure}

Using LATTICEEASY, the evolution of the Einstein-frame scalar field $\chi$ with time is obtained as shown in Fig. \ref{figure:chi-t-2}. We see in the figure that the $\chi$ field oscillates around the potential minimum which is located at $\chi = 0$, and as time evolves, the amplitude of the oscillations reduces. Having the time evolution of $\chi$ in hand, we are able to draw the diagram of the evolution of $M_{\Psi}^{2}$ and $\dot{\Omega}^{2}$ vs. time as illustrated in Fig. \ref{figure:MPsi2_Omegadot2-2}. It is evident from this figure that during the examined time interval, the inequality $\dot{\Omega}^{2}\ll M_{\Psi}^{2}$ is indeed valid, and therefore approximation \eqref{MPsit2-approx} that we used for $M_{\tilde{\Psi}}^{2}$ in our calculations, is completely viable.

\begin{figure}
\centering
\includegraphics[scale=0.4]{chi-t-2.pdf}
\caption{The variation of the Einstein-frame inflaton field ($\chi$) versus time in our setup.}
\label{figure:chi-t-2}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.4]{MPsi2_Omegadot2-2.pdf}
\caption{The evolution of $M_{\Psi}^{2}$ (black) and $\dot{\Omega}^{2}$ (orange) with time in our model.}
\label{figure:MPsi2_Omegadot2-2}
\end{figure}

The result of the LATTICEEASY code for the dispersion of the Einstein-frame preheat field, $\left\langle \tilde{\Psi}_{i}^2\right\rangle $, in terms of scale factor was presented in Fig. \ref{figure:avPsiti2-2}. Throughout this work, we normalize the scale factor to its value at the start of our simulation which is very close to the scale factor of the universe at the end of inflation. The figure shows that after a period of time, $\left\langle \tilde{\Psi}_{i}\right\rangle $ begins to grow an then reach its maximum value $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle \approx4\times10^{-5}$ at the scale factor $a\approx3.6$. After that, it decreases with scale factor, and for $a\gtrsim5$, it declines monotonically like $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle \propto a^{-2}$.

\begin{figure}
\centering
\includegraphics[scale=0.4]{avPsiti2-2.pdf}
\caption{The result of the LATTICEEASY code for the variation of the dispersion of the Einstein-frame preheat field, $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle $, in terms of scale factor in our framework.}
\label{figure:avPsiti2-2}
\end{figure}

By using the results for the time variation of $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle$ in Eq. \eqref{DeltaU}, we can calculate the potential correction term $\Delta U$ and include it in the background dynamics at each instant of time. Using the result for $\Delta U$ in Eq. \eqref{Utot}, we can evaluate the modified potential $U_{\mathrm{tot}}$. The time variation of $\Delta U$ in comparison with the one of the original potential $U$ has been depicted in Fig. \ref{figure:U_deltaU-t-2}. The figure implies that around the time $t\approx4.3\times10^{6}M_{P}^{-1}$, the correction term $\Delta U$ dominates over the original potential $U$, and its domination continues in all of the subsequent times. Due to this domination, the universe will undergo a transient accelerated expansion during preheating, before it settles to a decelerating phase again. The diagram of $\Delta U$ represents a peak with the value $\Delta U=1.277\times10^{-5}M_{P}^{4}$ at the time $t=4.915\times10^{6}M_{P}^{-1}$. At this moment, the vev of the preheat field squared gets the value of $\left\langle \tilde{\Psi}_{i}^{2}\right\rangle =3.021\times10^{-5}M_{P}^{2}$. By using this value in Eqs. \eqref{DeltaU} and \eqref{Utot}, we can draw the diagram of the potential correction term $\Delta U$ and the modified modified potential $U_{\mathrm{tot}}(\chi)$ at that moment as presented in Fig. \ref{figure:U-2}. In the figure, we have also compared the modified potential with the original potential $U(\chi)$. The effect of the correction term can be seen more clearly in the right panel of the figure. The figure shows obviously that the correction term $\Delta U$ lifts the potential so that it takes a positive value at its minimum while the original potential vanishes at that point.

\begin{figure}
\centering
\includegraphics[scale=0.4]{U_deltaU-t-2}
\caption{The evolution of the potential correction term $\Delta U$ (orange) and the original potential $U$ (black) versus time in the investigated model.}
\label{figure:U_deltaU-t-2}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.34]{U-2.pdf}
\caption{The modified potential $U_{\mathrm{tot}}(\chi)$ (solid line) in comparison with the original potential $U(\chi)$ (dash-dotted line). The correction term $\Delta U$ also has been shown in the figure by the dashed line. The plot in the left panel has been drawn in the usual linear scale, while the one in the right panel has been drawn in the logarithmic scale for a better demonstration.}
\label{figure:U-2}
\end{figure}

During a period of time around the instance $t=4.915\times10^{6}M_{P}^{-1}$ that the correction term takes rather considerable values, the effect of the lift in the inflaton potential becomes decisive in the background dynamics of the universe. In this period, the non-vanishing value at the potential minimum causes the universe to exit temporarily from the deceleration phase and enter an acceleration phase of expansion. The change of behavior from the deceleration to acceleration can be seen in Fig. \ref{figure:a-2} which demonstrates the variation of the scale factor with time in this case. The transition from deceleration to a temporary acceleration can also be seen more clearly in Fig. \ref{figure:q-2} that displays the plot of the deceleration parameter $q\equiv-a\ddot{a}/\dot{a}^{2}$ versus time. To avoid the oscillations of this parameter, here we have calculated the average of this parameter on each cycle of the $\chi$ field oscillation. The solid curve in the figure corresponds to the case with the inclusion of $\Delta U$, while the dashed curve is related to the case without the inclusion of the correction term. In the case that $\Delta U$ is not included in the evolution of the background, the deceleration parameter is always positive which indicates the expansion of the universe remains always decelerating. However, in the case that $\Delta U$ was included, the parameter $q$ can get negative values in an interval of time. Before the time $4.6\times10^{6}M_{P}^{-1}$, we have $q\approx0.5>0$ and the universe is in a phase of deceleration. Between the times $4.6\times10^{6}M_{P}^{-1}$ and $7.8\times10^{6}M_{P}^{-1}$, we have $q<0$ which indicates that the universe experiences an accelerating phase. After $t\approx7.8\times10^{6}M_{P}^{-1}$, the universe transits again to a deceleration phase with $q>0$. The occurrence of a phase of temporary acceleration after inflation may have notable consequences on the theoretical and observational implications. For example, this stage of nonthermal acceleration during preheating may be useful in the dilution of moduli and gravitino which otherwise lead to problems \cite{Felder:2000sf}.

\begin{figure}
\centering
\includegraphics[scale=0.4]{a-2.pdf}
\caption{The evolution of scale factor with time in our scenario. The solid line shows the result by including $\Delta U$, while the dashed line shows the one without including it.}
\label{figure:a-2}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.4]{q-2.pdf}
\caption{The deceleration parameter $q\equiv-a\ddot{a}/\dot{a}^{2}$ with respect to cosmic time in our model.}
\label{figure:q-2}
\end{figure}

In Fig. \ref{figure:nk-2}, we have examined the impact of the correction term to the inflaton potential on the number density \eqref{nk} of the particles generated during preheating. In this figure, the dashed and solid plots correspond to the case with the inclusion of $\Delta U$ and the case without the inclusion of it, respectively. If we compare these plots with each other, we see that in the presence of $\Delta U$, the number density takes smaller values, but the efficiency of preheating is still enough to explain the process of particle production after inflation. Therefore, we conclude that consideration of the potential correction term decreases the efficiency of the preheating mechanism to some extent.

\begin{figure}
\centering
\includegraphics[scale=0.4]{nk-2.pdf}
\caption{The variation of number density with time in our investigation. The dashed plot shows the number density without the inclusion of the corrected term $\Delta U$, while its contribution has already been involved in the solid plot.}
\label{figure:nk-2}
\end{figure}

At the end of this section, we examine the GWs generated during preheating in our scenario. The spectrum of GWs provided by LATTICEEASY for this model has been demonstrated in Fig. \ref{figure:OmegaGWh2-2}. The solid curve in this figure shows the spectrum for the case where $\Delta U$ is included, whereas the dashed one shows the spectrum without it. In the case where $\Delta U$ is included, the peak of the spectrum is of order $\Omega_{\mathrm{GW}}h^{2}\sim2.6\times10^{-12}$, but in the case where correction term was not included, the peak is of order $\Omega_{\mathrm{GW}}h^{2}\sim3.2\times10^{-11}$. So the amplitude of GWs produced during preheating is reduced by an order of magnitude. For both cases, the spectrum appears in the frequency range $f\sim6\times10^{8}-4\times10^{10}\,\mathrm{Hz}$. The range of the frequency and amplitude of the spectrum in the case without the correction term is almost similar to those values in the benchmark setup of preheating \cite{Easther:2006vd}, but the amplitude of the spectrum in the case with the involving it is one order of magnitude less than the result of the benchmark scenario. We see that the inclusion of the correction term reduces the amplitude of the stochastic GWs generated during preheating by one order of magnitude compared to the case without it. This reduction also originates from the temporary acceleration phase induced by the
\begin{figure}
\centering
\includegraphics[scale=0.4]{OmegaGWh2-2.pdf}
\caption{The LATTICEEASY result for the spectrum of the GWs produced during preheating in our setup. The solid plot shows the spectrum of the case with the inclusion of the potential term $\Delta U$, while the dashed plot corresponds to the case without involving it in the background dynamics. The spectrum in both cases has been evaluated at the time $t=10^{8}M_{P}^{-1}$.}
\label{figure:OmegaGWh2-2}
\end{figure}

\section{Conclusions}
\label{section:conclusions}

Light preheat fields coupled to the inflaton, can drain the energy of the inflaton during the deceleration phase. However as the inflaton energy is pumped to the preheat field and it obtains nonzero number densities, nonzero dispersion for the preheat field is developed that can affect the preheating phase. Here we investigated this issue in the context of non-symmetry-breaking non-minimal inflation. The M-flation scenario has motivations from the string theory, and it can resolve some of the problems of the chaotic inflationary models. However, despite its successes, M-flation suffers from several shortcomings that challenge its foundation. One way to remedy the problems of M-flation is the setup of non-$\mathbb{M}$-flation in which the inflaton field couples non-minimally to gravity. In this framework, we can improve the consistency of the inflationary observables with the recent CMB observations, in comparison with the minimal M-flation. The number of the D3-branes in this framework could be considerably less than the required number in the minimal setup, which ameliorates the issue of backreaction.

The non-$\mathbb{M}$-flation model, like the minimal M-flation scenario, has this interesting feature: it possesses scalar and gauge spectator fields. These fields are frozen classically during inflation, but they can play an important role in the dynamics of the universe in the preheating era after inflation.  The quantum vev of the spectator field squared leads to additional energy which lifts the potential energy of the inflaton field. We examined the effects of the modified potential on the dynamics of the post-inflationary universe in the setting of non-$\mathbb{M}$-flation. In this framework, we focused on the non-symmetry-breaking model in which the potential has only one minimum in its shape in both the Jordan and Einstein frames, and the potential is completely symmetric around this minimum in both frames. This model satisfies the current CMB constraints from the Planck 2018 measurements on the inflationary observables. In addition, since the number density can grow substantially in this model, therefore it provides an appealing preheating mechanism.

We used the LATTICEEASY code to compute the dispersion of the Einstein-frame preheat field, $\left\langle \tilde{\Psi}_{i}^2\right\rangle $, at each instant of time. The results of this quantity were used to evaluate the potential correction term $\Delta U$ and also the modified potential $U_{\mathrm{tot}}$. We included the contribution of $\Delta U$ in the inflaton potential and examined their effects on the preheating process. Our findings imply that around a special time during preheating, the correction term $\Delta U$ dominates over the original potential $U$, and it continues to be dominant in all of the subsequent times. Due to the correction term $\Delta U$, the inflaton potential is lifted so that it acquires a positive value at its minimum while the original potential vanishes at that point. The non-vanishing value at the potential minimum causes the universe to transit to a temporary acceleration phase of expansion during its evolution in the preheating era. This transient period of acceleration is a remarkable consequence coming from consideration of the quantum correction term in the post-inflationary dynamics. The occurrence of this phase of acceleration after inflation may have notable consequences. For example, this period of nonthermal acceleration during preheating may help in solving the moduli and gravitino problems \cite{Felder:2000sf}.  We furthermore examined the impact of the correction term to the inflaton potential on the number density of the particles generated during preheating in our setup. We showed that in the presence of $\Delta U$, the number density takes smaller values compared to the case without it, but the efficiency of preheating is still large enough to explain the process of particle production after inflation. Therefore, we conclude that consideration of the potential correction term reduces the efficiency of the preheating probes to some extent.

Finally, we examined the effect of the potential correction term on the stochastic GWs generated during preheating in our model. For this purpose, we applied LATTICEEASY to evaluate the spectrum of these GWs in our setup. Our results demonstrate that if the contribution of the potential correction term is included in the universe dynamics, then the peak of the GWs spectrum becomes of order $\Omega_{\mathrm{GW}}h^{2}\sim2.6\times10^{-12}$, while in the case without involving it, the peak is of order $\Omega_{\mathrm{GW}}h^{2}\sim3.2\times10^{-11}$. For both cases, the spectrum appears in the frequency range $f\sim6\times10^{8}-4\times10^{10}\,\mathrm{Hz}$. The range of the frequency and amplitude of the spectrum of the case without the correction term is almost similar to those values in the benchmark setup of preheating \cite{Easther:2006vd}, but the amplitude of spectrum in the case with the involving $\Delta U$ is one order of magnitude less than the one in the benchmark scenario. Therefore, we conclude that the inclusion of the correction term reduces the amplitude of the stochastic GWs generated during preheating one order of magnitude compared to the case without it.