# Evolution of the spin, spectrum and super-orbital period of the
ultraluminous X-ray pulsar M51 ULX7.

Murray Brightman  Cahill Center for Astrophysics, California Institute of
Technology, 1216 East California Boulevard, Pasadena, CA 91125, USA  Matteo
Bachetti  INAF-Osservatorio Astronomico di Cagliari, via della Scienza 5,
I-09047 Selargius (CA), Italy  Hannah Earnshaw  Cahill Center for
Astrophysics, California Institute of Technology, 1216 East California
Boulevard, Pasadena, CA 91125, USA  Felix Fürst  Quasar Science Resources SL
for ESA, European Space Astronomy Centre (ESAC), Science Operations
Departement, 28692 Villanueva de la Cañada, Madrid, Spain  Marianne Heida
European Southern Observatory, Garching, Germany  Gian Luca Israel  INAF-
Osservatorio Astronomico di Roma, via Frascati 33, 00078 Monteporzio Catone,
Italy  Sean Pike  Cahill Center for Astrophysics, California Institute of
Technology, 1216 East California Boulevard, Pasadena, CA 91125, USA  Daniel
Stern  Jet Propulsion Laboratory, California Institute of Technology,
Pasadena, CA 91109, USA  Dominic J Walton  Institute of Astronomy, Madingley
Road, Cambridge CB3 0HA, UK<EMAIL_ADDRESS>
###### Abstract

M51 ULX7 is among a small group of known ultraluminous X-ray pulsars (ULXP).
The neutron star powering the source has a spin period of 2.8 s, orbits its
companion star with a period of 2 days, and a super-orbital period of 38 days
is evident in its X-ray lightcurve. Here we present NuSTAR and XMM-Newton data
on the source from 2019 obtained when the source was near its peak brightness.
We detect the pulsations, having spun up at a rate of 3$\pm 0.5\times
10^{-10}$ s s-1 since they were previously detected in 2018. The data also
provide the first high-quality broadband spectrum of the source. We find it to
be very similar to that of other ULXPs, with two disk-like components, and a
high energy tail. When combined with XMM-Newton data obtained in 2018, we
explore the evolution of the spectral components with super-orbital phase,
finding that the luminosity of the hotter component drives the super-orbital
flux modulation. The inclination the disk components appear to change with
phase, which may support the idea that these super-orbital periods are caused
by disk precession. We also reexamine the super-orbital period with 3 years of
Swift/XRT monitoring, finding that the period is variable, increasing from
38.2$\pm 0.5$ days in 2018–2019 to 44.2$\pm 0.9$ days in 2020–2021, which
rules out alternative explanations for the super-orbital period.

††facilities: Swift(XRT), NuSTAR, XMM-Newton††software: NuSTARDAS, XMMSAS
(Gabriel et al., 2004), XSPEC (Arnaud, 1996), hendrics (Bachetti, 2015)

##  1 Introduction

M51 ULX7 was first detected as an X-ray source in the galaxies of M51 (NGC
5194/5) by the Einstein X-ray observatory (Palumbo et al., 1985). The source
was observed with an X-ray luminosity of $\sim 10^{39}$ erg s-1, which the
authors noted was brighter than any X-ray source in our own Galaxy, or M31,
and that it exceeded the Eddington luminosity of a 1 $M_{\odot}$ object. They
suggested this X-ray source, and others like it, could be powered by a neutron
star experiencing super-Eddington accretion in a non-spherically symmetric
accretion flow geometry or could indicate the presence of a more massive black
hole. These luminous X-ray sources would then come to be known as
ultraluminous X-ray sources, or ULXs (Fabbiano, 1989, see also recent reviews
by Kaaret et al. (2017) and Fabrika et al. (2021)). Thanks to the sensitivity
and spatial resolution of XMM-Newton and Chandra, we now know hundreds of ULXs
(e.g. Liu & Mirabel, 2005; Swartz et al., 2011; Walton et al., 2011; Earnshaw
et al., 2019; Kovlakas et al., 2020).

While the assumption that ULXs were powered by black holes gained the most
traction since their discovery by Einstein, the neutron star hypothesis was
eventually proven, at least for a handful of sources. This began in 2014 with
the detection of coherent pulsations from the ULX M82 X-2 by the NuSTAR
observatory (Bachetti et al., 2014), determining the source to be powered by a
neutron star, since black holes are incapable of producing such signals. This
was followed by NGC 5907 ULX (Israel et al., 2017a), NGC 7793 P13 (Israel et
al., 2017b; Fürst et al., 2017), NGC 300 ULX (Carpano et al., 2018), NGC 1313
X-2 (Sathyaprakash et al., 2019), and the subject of this paper, M51 ULX7
(Rodríguez Castillo et al., 2020). SMC X-3 (Tsygankov et al., 2017), Swift
J0243.6+6124 (Wilson-Hodge et al., 2018) and RX J0209.6-7427 (Vasilopoulos et
al., 2020b) also briefly became ULX pulsars, and another candidate pulsating
ULX was also recently reported in NGC 7793 (Quintin et al., 2021).

From timing analysis of XMM-Newton data, Rodríguez Castillo et al. (2020)
determined the spin period of M51 ULX7 to be $2.8$ s, and that the neutron
star was in a 2-day orbit with a $>8$ $M_{\odot}$ companion star, making it a
high-mass X-ray binary system (HMXB). The long-term, secular spin up rate was
also found to be $\sim 10^{-9}$ s s-1. Hu et al. (2021) and Vasilopoulos et
al. (2021) also found evidence for periodic dips in the Chandra X-ray light
curve that are associated with the 2-day binary orbital period which they
interpret as eclipses, implying that the orbit of the neutron star and its
donor star is seen at high inclination.

The spin period is similar to the $\sim 1$ s spin periods of M82 X-2, NGC 5907
ULX, and NGC 7793 P13. All four sources also share the common characteristic
of having periodic flux modulations in their long-term X-ray lightcurves, with
the latter over a range of 60–80 days (Walton et al., 2016; Fürst et al.,
2018; Brightman et al., 2019). The flux modulation from M51 ULX7 was found to
have a period of 38 days from Swift/XRT monitoring (Vasilopoulos et al.,
2020a; Brightman et al., 2020b). These periodic flux modulations, which in
most cases are longer than the orbital period of the system, have been
interpreted as precession of a large scale height disk, possibly caused by the
Lense-Thirring effect (e.g. Middleton et al., 2018). However, the origin of
these super-orbital periods is still a matter of debate.

While the aforementioned ULXs have now been determined to be powered by a
neutron star, it is still unknown what fraction, if any, are powered by black
holes. From an X-ray spectral standpoint, the neutron star powered ULXs appear
very similar to ones with unknown accretors, albeit among the hardest,
implying that the vast majority could be powered by neutron stars (Pintore et
al., 2017; Koliopanos et al., 2017; Walton et al., 2018a; Gúrpide et al.,
2021). The spectral shape consists of two disk-like components, a cooler one
which may come from the outer regions of an accretion disk or the photosphere
of an outflow (Qiu & Feng, 2021), and a hotter component, which may originate
from the inner regions of the accretion disk (Walton et al., 2018a), an
accretion curtain (Mushtukov et al., 2017), or Compton up-scattering
(Titarchuk, 1994). A high energy tail is also seen when NuSTAR data are
available, and appear to be associated with the pulsed component (Walton et
al., 2018a).

While NuSTAR has observed M51 ULX7 on two previous occasions with lower energy
coverage, the exposure time was either too short for a good quality spectrum
(Earnshaw et al., 2016), or the source was caught at a low flux (Brightman et
al., 2018). In NuSTAR Cycle 5, we obtained joint NuSTAR and XMM-Newton
observations of M51 ULX7, timed to occur at the peak of the periodic flux
modulation, with the aims of obtaining a high-quality broadband X-ray spectrum
of ULX7, modeling its emission components, and tracking its pulsations. We
present the results from these observations in this paper. We assume a
distance of 8.58$\pm$0.10 Mpc to M51, derived from the tip of the red giant
branch method (McQuinn et al., 2016).

##  2 X-ray data reduction

###  2.1 NuSTAR

NuSTAR (Harrison et al., 2013) observed M51 from 2019-07-10 05:56:09 to
2019-07-14 01:26:09 with an exposure of 169 ks (obsID 60501023002). We used
heasoft v6.28, nustardas v2.0.0 and caldb v20201101 to analyze the data. We
produced cleaned and calibrated events files using nupipeline with the
settings saacalc=3 saamode=OPTIMIZED tentacle=yes to account for enhanced
background during passages of the South Atlantic Anomaly (SAA), which reduced
the exposure time to 162 ks. We used nuproducts to produce spectral data,
including source and background spectra, and response files. A circular region
with a radius of 30″ was used to extract the source spectra. A circular region
with a radius of 100″ was used to extract the background spectra, taking care
to use extract the background from the same chip as the source. For timing
analyses, we used the heasoft tool barycorr to apply a barycentric correction
to the event times of arrival.

###  2.2 XMM-Newton

XMM-Newton (Jansen et al., 2001) observed M51 from 2019-07-11 10:47:26 to
2019-07-12 08:09:24.000 with an exposure of 77 ks (obsID 0852030101). We used
xmmsas v18.0.0 to analyze the data. We first identify periods of high
background by creating a lightcurve of the events in the 10–12 keV band,
creating good time intervals where the rate was less than 0.1 counts s-1 in
this band, leaving 69 ks of data. Events were selected with PATTERN$\leq 4$
for the pn and PATTERN$\leq 12$ for the MOS. A circular region with a radius
of 30″ was used to extract the source spectrum. A circular region with a
radius of 60″ was used to extract the background spectra, on the same chip as
the source and also in the galaxy in order to account for the soft diffuse
emission the source is embedded in (e.g. Earnshaw et al., 2016). Data from the
pn and both MOS instruments were extracted in this way. For timing analyses,
we used the xmmsas tool barycen to apply a barycentric correction to the event
times of arrival.

###  2.3 Swift

We used the online tool provided by the University of
Leicester111https://www.swift.ac.uk/user_objects/ (Evans et al., 2007, 2009)
to extract the Swift/XRT (Burrows et al., 2005) lightcurve of ULX7. All
products from this tool are fully calibrated and corrected for effects such as
pile-up and the bad columns on the CCD. We selected observations with target
IDs 11417, 30083 and 32017, and binned the lightcurve in time, with the
maximum bin size of 500 ks (5.79 days) and a minimum detection of 2.5$\sigma$.
The lightcurve is plotted in Figure 1, with the time of the 2019 NuSTAR and
XMM-Newton observation marked, as well as the 2018 XMM-Newton observations.

Figure 1: Swift/XRT lightcurve of ULX7 over the period May 2018 – May 2021
(black data points). Downward pointing arrows show 2.5$\sigma$ upper limits.
XMM-Newton and Chandra observation times are marked with red and blue lines
respectively.

##  3 Pulsation Analysis

We searched for pulsations using the fast $Z^{2}_{n}$ search implemented in
hendrics tool HENzsearch (Bachetti, 2015). This tool folds the data along a
grid of frequencies and frequency derivatives and calculates the $Z^{2}_{n}$
statistic starting from the folded profiles (see Huppenkothen et al., 2019;
Bachetti et al., 2021). We used $n=1$ (Rayleigh test) and folded profiles of
16 bins (adequate for using the $Z^{2}_{1}$ search with the binned
approximation).

We found a clear peak in the $f$-$\dot{f}$ plane in the XMM-Newton data,
however, no such signal was found in the lower-count NuSTAR data. We noted the
best solution from HENzsearch and used HENphaseogram to calculate the pulse
time-of-arrival (TOA) in 16 intervals during the observation. Following this
the graphical tool pintk in pint (Luo et al., 2021) was used to fit the TOAs
and get a timing solution and its uncertainty. We find $f=0.358842(5)$ Hz
($P=2.78674$ s), and $\dot{f}=-1.19(6)\times 10^{-8}$. We used the solution
from pint to phase tag the events using HENphasetag.

We did not find evidence for the orbit of the neutron star and its companion
when searching for a second derivative in the pulse frequency. This is likely
due to too few counts and only a fraction of the orbit being covered by the
XMM-Newton observation. Therefore we could not make a correction for the orbit
in the determination of the spin period. This means that the observed first
derivative in the pulse frequency is the sum of both the secular spin up of
the neutron star and the orbital motion. Rodríguez Castillo et al. (2020)
calculated the maximum delay/advance introduced when not correcting for the
orbital motion, which was of the order of 1 ms. In the following, when
inferring the secular spin up of the pulsar, we assume this 1 ms value to be
the absolute uncertainty of the period.

Using the above timing solution and the phase-tagged events files, we create
pulse profiles by binning the lightcurve in 16 equally sized phase bins, and
energy bins of 0.5–0.1 keV, 1-2 keV and 2–10 keV. We plot these in Figure 2.
Using these, we also calculate the pulse fraction in each bin, defined as the
amplitude of the pulse divided by the mean count rate. These are plotted in
Figure 3. Also shown are the pulse fractions as determined by Rodríguez
Castillo et al. (2020) from the 2018 XMM-Newton observing campaign. We derive
an upper limit of 74% to the pulse fraction in the 10–20 keV band from the
NuSTAR data by searching 170 Fourier frequencies, and assuming $\dot{P}=\pm
2\times\dot{P}_{max}$.

Figure 2: The pulse profile of ULX7 in energy bands of 0.5–10 keV (black),
0.5–0.1 keV (blue), 1-2 keV (green) and 2–10 keV (red). Figure 3: The pulsed
fraction of ULX7 versus energy from the 2019 XMM-Newton data (black), and in
comparison to the 2018 data from Rodríguez Castillo et al. (2020) (green, red,
and blue for each obsID, slightly offset in energy for clarity). The pulsed
fraction is highly variable and generally increases with energy.

##  4 Spectral analysis

We begin analyzing the 2019 NuSTAR Cycle 5 NuSTAR and XMM-Newton spectral data
on ULX7 by grouping the spectra with a minimum of 1 count per bin. We load the
FPMA, FPMB, pn, MOS1 and MOS2 source spectra into xspec, subtracting the
background, and consider energies 3–20 keV for the NuSTAR data (the background
dominates the source above 20 keV so we do not consider these data), and
0.2–10 keV for the XMM-Newton data. We use the C-statistic, suitable for the
low number of counts per bin here and we use a constant term to take account
of cross-calibration uncertainties between instruments, which are typical of
those found in (Madsen et al., 2015).

Rodríguez Castillo et al. (2020) carried out fits to the series of high-
quality XMM-Newton spectra they obtained in 2018 when the pulsations were
detected. Their best-fit model consisted of an absorbed disk black body model
plus a hotter blackbody model (tbabs*ztbabs*(diskbb+bbodyrad)), which we fit
here. The tbabs model accounts for absorption in our Galaxy, fixed at
$3.19\times 10^{20}$ cm-2 (HI4PI Collaboration et al., 2016) and the ztbabs
accounts for absorption at the redshift of M51, $z=0.002$, left as a free
parameter. This model resulted in $C=4027.81$ with 4286 d.o.f.

However, there appears to be an excess at energies above 10 keV when fitting
with this model to the 2019 NuSTAR+XMM-Newton data (Figure 4). Rodríguez
Castillo et al. (2020) noted that in Walton et al. (2018a) a third component,
a powerlaw with a high energy cut off (cutoffpl) is used to model this
component which is attributable to the accretion column. This pulsed component
was isolated in the ULX pulsars M82 X-2, NGC 7793 P13 and NGC 5907 ULX by
extracting spectra from the brightest and the faintest quarters
($\Delta\phi_{\rm pulse}=0.25$) of the pulse cycle and subtracting the latter
from the former (i.e., pulse on $-$ pulse off, Brightman et al., 2016; Walton
et al., 2018b, a).

Figure 4: Top - XMM-Newton (black, red, and green), and NuSTAR (blue) spectra
of M51 ULX7. The spectra can be modeled by the typical ULX spectrum which
consists of two disk-like components, one at lower energy, the other at higher
energy. A third component seen at the highest energies is the pulsed
component, isolated using phase-resolved spectroscopy, defined in xspec as
tbabs*(diskbb+diskpbb+cutoffpl). Bottom - Data to model residuals for the
models described in Section 4.

We do the same pulse on $-$ pulse off analysis here for the XMM-Newton data on
ULX7 using the timing solution of the pulses found in Section 3, and fit it
with tbabs*ztbabs*cutoffpl, where the absorption components have been fixed to
the values from the time averaged spectrum. We cannot constrain $E_{\rm cut}$
due to the lack of high energy coverage, so we fix $E_{\rm cut}=8.1$, the
average values from Walton et al. (2018a). We can constrain $\Gamma=0.8\pm
0.3$ with $C=881.89$ with 946 d.o.f., which is consistent with the average
value of 0.5 from Walton et al. (2018a). We plot this spectrum in Figure 5.

Figure 5: Pulsed spectrum of ULX7 produced by the pulse on $-$ pulse off
analysis of the XMM-Newton data.

We then add this component to the time averaged spectral fit, keeping $\Gamma$
and $E_{\rm cut}$ fixed, but leaving the normalization free. This results in
$C=3954.59$ with 4285 d.o.f., presenting an improvement of the overall time-
averaged fit.

We note that if a diskpbb model were used for the hottest disk component
instead of bbodyrad, consistent with the models used in Walton et al. (2018a),
we get a better fit of $C=3951.17$ with 4284 d.o.f. Therefore, we present this
diskbb+diskpbb+cutoffpl model as our best fit.

Table 1: Results from the X-ray spectral modeling of joint XMM-Newton and NuSTAR data on ULX7 Parameter | Model
---|---
| diskbb+bbodyrad | diskbb+bbodyrad+cutoffpl | diskbb+diskpbb+cutoffpl
$N_{\rm{H}}$/$10^{20}$cm-2 | 1.3$\pm$0.7 | 3.3$\pm$0.7 | 3.9${}^{+1.3}_{-0.6}$
$T_{\rm in,cool}$/keV | 0.46$\pm 0.02$ | 0.36$\pm 0.02$ | 0.33$\pm 0.03$
$N_{\rm cool}/10^{-4}$ a | 4050${}^{+500}_{-480}$ | 950${}^{+350}_{-250}$ | 13500${}^{+6200}_{-1800}$
$T_{\rm hot}$/keV | 1.64$\pm 0.05$ | 1.39$\pm 0.09$ | 2.26${}^{+0.54}_{-0.17}$
$p$ | - | - | 0.99${}^{+u}_{-0.25}$
$N_{\rm hot}/10^{-4}$ a | 75.6${}^{+10.6}_{-9.5}$ | 71.2${}^{+17}_{-15}$ | 17.2${}^{+17}_{-12}$
$F_{\rm{X}}$/$10^{-13}$ erg cm-2 s-1 b | 8.3$\pm 0.3$ | 8.2$\pm 0.3$ | 7.7$\pm 0.2$
$L_{\rm{X}}$/$10^{39}$ erg s-1 c | 7.7$\pm 0.4$ | 7.6$\pm 0.4$ | 7.8$\pm 0.4$
$C_{\rm pn}$ | 1.00${}^{+0.07}_{-0.06}$ | 0.99$\pm 0.07$ | 0.99${}^{+0.07}_{-0.06}$
$C_{\rm MOS1}$ | 1.03${}^{+0.08}_{-0.07}$ | 1.04${}^{+0.08}_{-0.07}$ | 1.03${}^{+0.08}_{-0.07}$
$C_{\rm MOS2}$ | 1.05${}^{+0.08}_{-0.07}$ | 1.07${}^{+0.08}_{-0.07}$ | 1.05${}^{+0.08}_{-0.07}$
$C_{\rm FPMB}$ | 1.00$\pm 0.08$ | 1.02${}^{+0.09}_{-0.08}$ | 1.00$\pm 0.08$
$C$-statistic/d.o.f. | 4027.81/4286 | 3954.59/4285 | 3951.17/4284

##  5 Evolution of the spectral parameters

With the newly determined best-fit model found in Section 4, we proceed to
apply it to the high-quality data from XMM-Newton in 2018 in order to explore
the evolution of the X-ray spectral parameters with luminosity and super-
orbital phase. While these data were presented by Rodríguez Castillo et al.
(2020), they used a slightly different model, so we cannot strictly compare
their results to ours. We note that these models are purely phenomenological,
and any physical interpretations should have this added caveat.

We reduce the XMM-Newton obsIDs 0824450901, 0830191401, 0830191501, 0830191601
in the same way as described in Section 2.2. We then fit the spectra with the
diskbb+diskpbb+cutoffpl model described above. As found by Rodríguez Castillo
et al. (2020), ULX7 was seen in a low flux state in obsID 0830191401, and only
the low temperature diskbb can be seen. We remove the cutoffpl model from the
fit since it clearly over fits the data, and we fix the temperature and $p$
parameter of the diskpbb model to typical values in order to place an upper
limit on the normalization and flux of this component.

For each epoch, we calculate the total intrinsic luminosity and the luminosity
of the two disk components separately, determined using the model component
cflux and calculated in the 0.3–10 keV range. To get the super-orbital phase,
we divide the number of days since 2018-05-01 by the period of 38.3 days, and
subtract off the nearest integer. We also adjust the phase so that phase=0
corresponds to the peak of the super-orbital modulation. In Brightman et al.
(2020b), we showed that the super-orbital period of 38.3 days was consistent
across the 2018–2019 epoch we study here.

First we plot the total intrinsic luminosity, and the luminosity of the two
disk components as a function of super-orbital phase, and compare these to the
average total luminosity as seen from the Swift/XRT monitoring. Figure 6 shows
that the total luminosity measured by XMM-Newton agrees in general with that
of the average luminosity seen by Swift/XRT. In terms of the two disk
components, the cool component’s luminosity remains relatively constant,
whereas the hot component’s luminosity varies more, albeit with large
uncertainties.

We also plot the pulsed fraction derived in Section 3 and in Rodríguez
Castillo et al. (2020), and how it varies with super-orbital phase in Figure
6. We see a similar picture here, where the pulsed fraction in the low energy
band varies relatively little over the super-orbital period, but the higher
energy pulsed fraction varies strongly. This is likely driven by the variation
in the non-pulsed, hot disk component.

Figure 6: Top - The total luminosity (black squares) and luminosity of the
cool (blue squares) and hot (red squares) disk components of ULX7 as a
function of its super-orbital phase from NuSTAR and XMM-Newton observations in
2018–2019 (offset in phase for clarity). The average luminosity as a function
of phase from Swift/XRT observations is shown with a black line. Bottom - The
pulsed fraction in the 0.5–1 keV band (blue), 1–2 keV band, and 2–10 keV band
(red) as a function of super-orbital phase with data from Rodríguez Castillo
et al. (2020).

We then proceed to plot the spectral parameters of the disk components against
the intrinsic luminosity and super-orbital phase in Figure 7. We do not see
any statistically significant correlations between the spectral parameters and
$L_{\rm{X}}$ for either the cool diskbb component or hot diskbb component. We
do see that $N_{\rm{H}}$, $T_{\rm in,cool}$, $N_{\rm cool}$ and $N_{\rm hot}$
appear to vary with the phase of the super-orbital period, however we likely
do not have enough data to claim a significant relationship. If the dependence
is real, then since the normalization of the disk components depends on the
orientation of the disk, this may imply that the inclination of the disks os
changing with super-orbital phase, or in other words, the disks are
precessing.

However, since the disk temperature and normalization parameters are
degenerate with each other, it is possible that these variations are driven by
the degeneracy instead. We investigate this possibility by using Markov-Chain
Monte-Carlo methods to map out the $T_{\rm in,cool}$ and $N_{\rm cool}$
parameter space. In xspec we use the Goodman-Weare algorithm with 8 walkers
and a total length of 10000 steps with a burn-in phase of 5000 steps. Figure 8
shows the results of this. While a degeneracy can be seen between the two
parameters, the results map out regions which are almost mutually exclusive,
implying that the degeneracy is not the cause of the possible dependence of
these parameters on super-orbital phase.

Figure 7: Spectral parameters of the diskbb (cool component, blue data points)
and diskpbb (hot component, red data points) models fitted to the 2019 NuSTAR
and XMM-Newton data, and 2018 XMM-Newton data, plotted against the
$L_{\rm{X}}$ for each individual component (left, 0.3–10 keV unabsorbed) and
super-orbital phase (right, two cycles plotted for clarity). Uncertainties are
plotted at the 90% confidence level. Figure 8: Results of a 10000-step MCMC
analysis for the $T_{\rm in,cool}$ and $N_{\rm cool}$ parameters from each
XMM-Newton observation. While a degeneracy can be seen between the two
parameters, the results map out regions which are almost mutually exclusive,
implying that the degeneracy is not the cause of the possible dependence of
these parameters on super-orbital phase.

##  6 Evolution of the super-orbital period

In Brightman et al. (2020b), we showed that ULX7 exhibited a super-orbital
period of 38 days over 500 days ($\sim$1.5 years) from May 2018 (see also
Vasilopoulos et al., 2020a). However, when we applied our analysis to the most
recent, full $\sim 3$-year data set, the signal appears much weaker, with a
peak in the $L$ statistic at 34.1 days and a peak in the Lomb-Scargle
periodogram at 39.0 days. This implies that after 500 days, the super-orbital
period disappears, or changes phase and or period. Figure 9 shows the full
3-year lightcurve with the average profile of the flux modulations from the
first 500 days overplotted. The deviation from the profile after 500 days is
clear in the residuals, and exhibits both positive and negative deviations
therefore cannot be explained by the anomalous low flux states seen by
Vasilopoulos et al. (2020a).

We investigate further by splitting the lightcurve into smaller sections in
time to determine how the super-orbital period evolves. As in Brightman et al.
(2020b), we use epoch folding on the $L$ statistic (Davies, 1990), using 10
phase bins. This time however, since we know the approximate period, we search
over a narrower range of 30–50 days in 200 equally spaced bins for the epoch
folding. We find that a 300-day section is sufficient to recover the 38-day
period at the beginning of the lightcurve, and then we progress through the
lightcurve with steps of 30 days.

In Brightman et al. (2020b), we used simulations to determine the false alarm
rate for the 38-day period signal, finding it to be $>99.9$% significant. We
repeat this analysis for our 300-day bins, simulating 10,000 lightcurves with
2000 s resolution and a red noise power spectrum and we sample them with the
same observational sampling as the real lightcurves. We then note the largest
peak in each periodogram, irrespective of period. We define the false alarm
rate as the number of simulated lightcurves that produce a peak as high as the
real one, divided by the total number of simulations. We do this for each of
the time bins since each time bin has a different number and spacing of
observations. For each time bin, we determine the $L$-stat level which
corresponds to a 0.3% false alarm rate, equivalent to a 3$\sigma$ detection.
In Figure 10 we plot the super-orbital periods which are detected at
$>3\sigma$ against time.

In order to determine the uncertainties on the periods detected above, we use
a Monte-Carlo technique to resample the lightcurve where the count rate of
each observation is a random number that is taken from a Gaussian distribution
with a mean equal to the observed count rate and a standard deviation equal to
the uncertainty on the count rate. We do this 10,000 times, rerunning our
analysis on each lightcurve, and noting the peak period. We define the
1$\sigma$ uncertainty as the standard deviation of periods recovered. We plot
these derived uncertainties in Figure 10.

We find that at the beginning of the lightcurve the super-orbital period is
detected at $>3\sigma$ at $\sim$38 days, as found in Brightman et al. (2020b),
but at the end of the lightcurve it is detected at $\sim$44 days. The highest
peaks in the periodograms are at 0–300 days, with $P=38.2\pm 0.5$ days, and
780–1080 days where $P=44.2\pm 0.9$ days. These time bins are independent and
the findings indicate that the super-orbital period has increased in length.
There are no obvious correlations between the super-orbital period and the
number of observations per bin, or the average count rate.

While the above approach accounts for the statistical uncertainties on the
individual measurements, it does not account for possible effects due to
sampling. In order to do this we employ a bootstrapping method, whereby we
randomly exclude 10% of the observations in each time bin. Repeating this
10,000 times, we find standard deviations on the periods derived were 0.13 for
the 38.2 day period 0.58 for the 44.2 day period. If we increase the number of
excluded pointings to 20%, this becomes 0.28 and 1.83 respectively. This
implies the statistical uncertainty on the measurements dominates the
uncertainty in the super-orbital period.

Finally, we also check if the longer super-orbital period seen at 780–1080
days can be produced by the 38.2-day signal, but with the lightcurve sampling
at 780–1080 days. We do this by projecting the 38.2 day profile to 780–1080
days and simulating a lightcurve with the same observing strategy during that
time. Again, we do this 10,000 times, however none of these simulations
produce a peak at 44.2-days or greater with a height at least equal to the
observed periodogram. It is therefore highly unlikely that the Swift observing
strategy caused the observed change in super-orbital period.

We show the lightcurves and folded lightcurves from the 0–300 days and
780–1080 days epochs in Figure 11. The profiles from the two lightcurves
appear similar, albeit perhaps with most recent one being more peaked, and
less sinusoidal.

Figure 9: Top - Swift/XRT lightcurve of ULX7 over 3-years 2018–2021 (black
data points) with the average profile of the 38.2-day super-orbital period
seen in the first 500 days overplotted (red line). Bottom - Residuals of the
data to the profile. The data clearly follow the average profile for the first
500 days, but deviate thereafter, which we find is in part due to a change of
the super-orbital period. Figure 10: Left - The evolution of the super-orbital
period, where each data point represents a 3$\sigma$ detection for a 300-day
time bin. The detection significance is determined by simulations. Error bars
are 1$\sigma$, and are also determined by simulations. The super-orbital flux
modulation is detected at $>3\sigma$ at the beginning and end of the
lightcurve, increasing in period from $\sim$38 days to $\sim$44 days. Middle -
The number of Swift/XRT observations used in each time bin. Right - The
average Swift/XRT count rate observed in each time bin.

Figure 11: Swift/XRT lightcurves and folded lightcurves over the 0–300 days and 780–1080 days epochs, showing the change in super-orbital period. Table 2: Details of the observations used in this work Observatory | ObsID | Start time | Exposurea | Count rate | Fluxb | Super-orbital phase
---|---|---|---|---|---|---
|  | (UT) | (ks) | (counts s-1) | ( erg cm-2 s-1) |
XMM-Newton | 0824450901 | 2018-05-13 21:18:47 | 60.6 | 0.167$\pm 0.002$ | 5.8$\times 10^{-13}$ | 0.15
XMM-Newton | 0830191401 | 2018-05-25 20:26:58 | 77.6 | 0.018$\pm 0.001$ | 2.9$\times 10^{-14}$ | -0.53
XMM-Newton | 0830191501 | 2018-06-13 01:39:03 | 48.2 | 0.221$\pm 0.002$ | 7.9$\times 10^{-13}$ | -0.04
XMM-Newton | 0830191601 | 2018-06-15 01:24:21 | 49.3 | 0.209$\pm 0.002$ | 7.5$\times 10^{-13}$ | 0.01
NuSTAR | 60501023002 | 2019-07-10 05:56:09 | 162 | 0.006$\pm 0.0002$ | 7.6$\times 10^{-13}$ | 0.26
XMM-Newton | 0852030101 | 2019-07-11 10:47:26 | 58.8 | 0.209$\pm 0.002$ | 7.6$\times 10^{-13}$ | 0.26

##  7 Discussion

###  7.1 The spin of the neutron star

The pulsations detected with $P=2.78674$ s in our 2019 XMM-Newton dataset of
M51 ULX7 represent a secular spin up of the neutron star from $P=2.79812$ s in
2018 (Rodríguez Castillo et al., 2020), with an average spin up rate of
3$\times 10^{-10}$ s s-1 over the 425-days spanning the observations. As
described in Section 3, we were not able to correct the 2019 data for the
orbital motion of the pulsar. We determined that this produced an uncertainty
of order 1 ms in the period. This then leads to an uncertainty of 0.5$\times
10^{-10}$ s s-1 in the derived secular spin-up rate. The secular spin up rate
is similar to that calculated for the short timescale covered by the 2018
data, and appears comparable with the spin up of $\sim 10^{-9}$ s s-1
calculated from observations spanning 13 years, which is in the range of other
ULX pulsars (Rodríguez Castillo et al., 2020) .

###  7.2 The broadband X-ray spectrum

We presented here the first high-quality broadband X-ray spectrum with
coverage above 10 keV of M51 ULX7, obtained with NuSTAR and XMM-Newton, over
the energy range 0.2–20 keV. Other members of the known ULX pulsar population
with at least one high-quality broadband X-ray spectrum are NGC 5907 ULX
(Walton et al., 2015), NGC 7793 P13 (Walton et al., 2018b), and NGC 300 ULX
(Carpano et al., 2018). M82 X-2 has broadband spectra but with significant
confusion due to the proximity of M82 X-1 (Brightman et al., 2020a), and NGC
1313 X-2 is not well detected above 10 keV (Bachetti et al., 2013). Therefore
adding to this small sample of ULX pulsars with a high-quality broadband X-ray
spectrum is significant. Indeed, the sample of ULXs with high quality
broadband spectra in general, be it with a known neutron star accretor or an
unknown accretor, is also small (Walton et al., 2018b).

The broadband X-ray spectrum of M51 ULX7 can be well described by two disk-
like components, plus a higher energy component associated with the
pulsations. This is qualitatively very similar to the other ULX pulsars, and
indeed, as pointed out in Walton et al. (2018a), very similar to all ULXs in
general, regardless of the known or unknown accretor (see also Pintore et al.,
2017; Koliopanos et al., 2017). As presented by Walton et al. (2018a), the
temperature of the cool component for ULXs ranges from 0.2–0.5 keV. We find
for ULX7 that it is 0.33 keV. The hotter component has a temperature range of
1.2–3 keV, and is 2.5 keV for ULX7, making ULX7 indistinguishable from other
ULXs in terms of its disk temperatures. Walton et al. (2018a) also noted that
the temperature ratio for ULX pulsars in their sample was $\sim 3$, while the
other ULXs had a temperature ratio of $\sim 8$. For the ULX7, the ratio is
7.5, therefore more in line with the other ULXs, rather than the other ULX
pulsars. Walton et al. (2018a) also found that the flux ratio of the pulsed
component, modeled by cutoffpl to the total flux in the 0.3–40 keV band was
higher for ULX pulsars than for the other ULXs. For M51 ULX7, however, this
ratio is 0.15, which is relatively low in comparison to NGC 5907 ULX (0.82 in
the high state), and NGC 7793 P13 (0.59).

###  7.3 The spectral evolution

The spectral evolution of ULX7 has been studied before. Both Yoshida et al.
(2010) and Earnshaw et al. (2016) fitted the available Chandra and XMM-Newton
data at the time with a power-law model. Earnshaw et al. (2016) found that the
spectral slope did not change significantly, despite the large changes in
$L_{\rm{X}}$. This included the very faint states observed with Chandra.
Rodríguez Castillo et al. (2020) also investigated the evolution of the
spectral parameters in their high-quality XMM-Newton data from 2018, using
their disk+black body model, also finding limited evidence for spectral
variations, as did Gúrpide et al. (2021).

In our investigation of the spectral parameters of ULX7 and how they depend on
the super-orbital phase of the system, we have found that the flux modulations
are primarily driven by changes in the flux of the hotter disk-like component,
with the flux of the cooler disk component changing less dramatically. We also
see a potential dependence of both the disk temperature and normalization of
the cooler disk-like component on super-orbital phase, even though the flux of
this component does not change significantly with the super-orbital phase.
However we do not have enough data to claim a significant relationship.

Since the normalization of the diskbb and diskpbb components are directly
related to the inclination of the disk (also the distance to the source and
inner disk radius), this potential dependence could be straightforwardly
interpreted as disk precession, which has been suggested as the mechanism for
the super-orbital periodic flux variations seen in ULX pulsars (e.g. Fürst et
al., 2017; Dauser et al., 2017; Middleton et al., 2018). Indeed, Fürst et al.
(2017) also found evidence that the spectral parameters of the disk-like
component of NGC 5907 ULX were dependent on the super-orbital phase of the
system. In this case it was the hot disk-like component, rather than the cool
one we see it in, and it was the radial temperature index, $p$, rather than
the normalization/temperature. Since NGC 5907 ULX is relatively absorbed, the
cool disk component is not visible. We note, however, that the models used are
purely phenomenological, and any physical interpretations should have this
added caveat.

The lack of any strong evolution of $N_{\rm{H}}$ with super-orbital period
rules out a warped accretion disk which periodically obscures the X-ray
source. While small variations in $N_{\rm{H}}$ can be seen, much larger
changes are needed to produce the $>1$ magnitude flux modulation.

###  7.4 The evolution of the super-orbital period

Having discovered that the X-ray flux from M51 ULX7 is modulated on a period
of 38.2$\pm 0.5$ days in Brightman et al. (2020b) from 2018–2019, we find here
that period is variable, and has shifted to 44.2$\pm 0.9$ days in 2020–2021.
Vasilopoulos et al. (2020a) also explored the Swift/XRT data on ULX7, finding
the same 38-day period and also found tentative evidence for a 49 day period
from Swift/XRT data taken in 2011, albeit covering only 1–2 cycles and
therefore quite uncertain.

An example of a super-orbital flux modulation in a source where the period is
variable is SMC X-1, which has a super-orbital period of $\sim 60$ days
(Gruber & Rothschild, 1984), and exhibits recurrent excursions to shorter
super-orbital periods, which may be (quasi)periodic themselves (Hu et al.,
2019). Unlike M51 ULX7, SMC X-1 shows strong spectral variability across its
super-orbital phase, with a harder spectrum at low fluxes, suggestive of
absorption being the main cause. This absorption appears weaker during the
super-orbital period excursion. From modeling of the pulse profile with super-
orbital phase from XMM-Newton and NuSTAR data, Brumback et al. (2020) found
that the pulse shape and phase of SMC X-1 are consistent with reprocessed
emission from a precessing inner disk. The wind-fed HMXB IGR J16493-4348 also
exhibits a variable super-orbital period, however in that case, the amplitude
is found to be variable, rather than the period. They suggest this is linked
to a variable accretion rate. However, they note that the timing and spectral
properties of ULXs show significant differences compared to those observed in
wind-fed HMXBs such as IGR J16493-4348.

As discussed above, the evolution of the spectral parameters across super-
orbital phase support the hypothesis that the flux modulations from M51 ULX7
are caused by precession of the accretion disk. Other theories such as a
warped accretion disk periodically obscuring the X-ray source as suggested for
SMC X-1 can be ruled out by the lack of $N_{\rm{H}}$ variations. It has also
been suggested that a third orbiting star could be causing the variations
(e.g. Middleton et al., 2018; Rodríguez Castillo et al., 2020). This now
appears unlikely since if the orbit of the third star were causing the
periodic variations, the orbit would have needed to change significantly in a
short time.

Vasilopoulos et al. (2020a) discussed the potential for free-precession of the
neutron star, which has been invoked to explain long term periodic changes of
isolated neutron stars. However, they note that free precession of the neutron
star alone cannot account for the variation in super-orbital period, which we
find here.

##  8 Conclusions

In our new NuSTAR and XMM-Newton data on M51 ULX7, we have found that the
neutron star powering the source has spun up at a rate of 3$\pm 0.5\times
10^{-10}$ s s-1 since the previous observations by XMM-Newton, which is
similar to that seen in other ULX pulsars. The data also provide the first
high-quality broadband spectrum, consisting of two disk-like components, the
temperatures of which are indistinguishable from other ULXs, and a high energy
tail. We found that the luminosity of the hotter component drives the super-
orbital flux modulation seen from the source. Finally we discovered that the
super-orbital period varies, and has increased from 38.2$\pm 0.5$ days in
2018–2019 to 44.2$\pm 0.9$ days in 2020–2021. This change in period rules out
some alternative explanations of the super-orbital period.

We thank the referee for their thorough review of our paper which improved it.
This work was also supported under NASA Contract No. NNG08FD60C. NuSTAR is a
project led by the California Institute of Technology, managed by the Jet
Propulsion Laboratory, and funded by the National Aeronautics and Space
Administration. This research has made use of the NuSTAR Data Analysis
Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC,
Italy) and the California Institute of Technology (USA). This work was also
based on observations obtained with XMM-Newton, an ESA science mission with
instruments and contributions directly funded by ESA Member States and NASA We
wish to thank the Swift PI, Brad Cenko for approving the target of opportunity
requests we made to observe M51, as well as the rest of the Swift team for
carrying them out. We also acknowledge the use of public data from the Swift
data archive. This work made use of data supplied by the UK Swift Science Data
Centre at the University of Leicester. This research has made use of data
and/or software provided by the High Energy Astrophysics Science Archive
Research Center (HEASARC), which is a service of the Astrophysics Science
Division at NASA/GSFC. GLI acknowledges funding from the Italian MIUR PRIN
grant 2017LJ39LM. The work of DS was carried out at the Jet Propulsion
Laboratory, California Institute of Technology, under a contract with NASA.
DJW acknowledges support from an STFC Ernest Rutherford Fellowship.

## References

  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17–+
  * Bachetti (2015) Bachetti, M. 2015, MaLTPyNT: Quick look timing analysis for NuSTAR data, ascl:1502.021
  * Bachetti et al. (2021) Bachetti, M., Pilia, M., Huppenkothen, D., et al. 2021, ApJ, 909, 33
  * Bachetti et al. (2013) Bachetti, M., Rana, V., Walton, D. J., et al. 2013, ApJ, 778, 163
  * Bachetti et al. (2014) Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202
  * Brightman et al. (2020a) Brightman, M., Walton, D. J., Xu, Y., et al. 2020a, ApJ, 889, 71
  * Brightman et al. (2016) Brightman, M., Harrison, F., Walton, D. J., et al. 2016, ApJ, 816, 60
  * Brightman et al. (2018) Brightman, M., Baloković, M., Koss, M., et al. 2018, ApJ, 867, 110
  * Brightman et al. (2019) Brightman, M., Harrison, F. A., Bachetti, M., et al. 2019, ApJ, 873, 115
  * Brightman et al. (2020b) Brightman, M., Earnshaw, H., Fürst, F., et al. 2020b, ApJ, 895, 127
  * Brumback et al. (2020) Brumback, M. C., Hickox, R. C., Fürst, F. S., et al. 2020, ApJ, 888, 125
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Science Reviews, 120, 165
  * Carpano et al. (2018) Carpano, S., Haberl, F., Maitra, C., & Vasilopoulos, G. 2018, MNRAS, 476, L45
  * Dauser et al. (2017) Dauser, T., Middleton, M., & Wilms, J. 2017, MNRAS, 466, 2236
  * Davies (1990) Davies, S. R. 1990, MNRAS, 244, 93
  * Earnshaw et al. (2016) Earnshaw, H. M., Roberts, T. P., Heil, L. M., et al. 2016, MNRAS, 456, 3840
  * Earnshaw et al. (2019) Earnshaw, H. P., Roberts, T. P., Middleton, M. J., Walton, D. J., & Mateos, S. 2019, MNRAS, 483, 5554
  * Evans et al. (2007) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379
  * Evans et al. (2009) —. 2009, MNRAS, 397, 1177
  * Fabbiano (1989) Fabbiano, G. 1989, ARA&A, 27, 87
  * Fabrika et al. (2021) Fabrika, S. N., Atapin, K. E., Vinokurov, A. S., & Sholukhova, O. N. 2021, Astrophysical Bulletin, 76, 6
  * Fürst et al. (2017) Fürst, F., Walton, D. J., Stern, D., et al. 2017, ApJ, 834, 77
  * Fürst et al. (2018) Fürst, F., Walton, D. J., Heida, M., et al. 2018, A&A, 616, A186
  * Gabriel et al. (2004) Gabriel, C., Denby, M., Fyfe, D. J., et al. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 314, Astronomical Data Analysis Software and Systems (ADASS) XIII, ed. F. Ochsenbein, M. G. Allen, & D. Egret, 759
  * Gruber & Rothschild (1984) Gruber, D. E., & Rothschild, R. E. 1984, ApJ, 283, 546
  * Gúrpide et al. (2021) Gúrpide, A., Godet, O., Koliopanos, F., Webb, N., & Olive, J. F. 2021, A&A, 649, A104
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., Hailey, C. J., & Zhang, W. W. 2013, ApJ, 770, 103
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hu et al. (2019) Hu, C.-P., Mihara, T., Sugizaki, M., Ueda, Y., & Enoto, T. 2019, ApJ, 885, 123
  * Hu et al. (2021) Hu, C.-P., Ueda, Y., & Enoto, T. 2021, ApJ, 909, 5
  * Huppenkothen et al. (2019) Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, ApJ, 881, 39
  * Israel et al. (2017a) Israel, G. L., Belfiore, A., Stella, L., et al. 2017a, Science, 355, 817
  * Israel et al. (2017b) Israel, G. L., Papitto, A., Esposito, P., et al. 2017b, MNRAS, 466, L48
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1
  * Kaaret et al. (2017) Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303
  * Koliopanos et al. (2017) Koliopanos, F., Vasilopoulos, G., Godet, O., et al. 2017, A&A, 608, A47
  * Kovlakas et al. (2020) Kovlakas, K., Zezas, A., Andrews, J. J., et al. 2020, MNRAS, 498, 4790
  * Liu & Mirabel (2005) Liu, Q. Z., & Mirabel, I. F. 2005, A&A, 429, 1125
  * Luo et al. (2021) Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45
  * Madsen et al. (2015) Madsen, K. K., Harrison, F. A., Markwardt, C. B., et al. 2015, ApJS, 220, 8
  * McQuinn et al. (2016) McQuinn, K. B. W., Skillman, E. D., Dolphin, A. E., Berg, D., & Kennicutt, R. 2016, ApJ, 826, 21
  * Middleton et al. (2018) Middleton, M. J., Fragile, P. C., Bachetti, M., et al. 2018, MNRAS, 475, 154
  * Mushtukov et al. (2017) Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Ingram, A. 2017, MNRAS, 467, 1202
  * Palumbo et al. (1985) Palumbo, G. G. C., Fabbiano, G., Trinchieri, G., & Fransson, C. 1985, ApJ, 298, 259
  * Pintore et al. (2017) Pintore, F., Zampieri, L., Stella, L., et al. 2017, ApJ, 836, 113
  * Qiu & Feng (2021) Qiu, Y., & Feng, H. 2021, ApJ, 906, 36
  * Quintin et al. (2021) Quintin, E., Webb, N. A., Gúrpide, A., Bachetti, M., & Fürst, F. 2021, MNRAS, 503, 5485
  * Rodríguez Castillo et al. (2020) Rodríguez Castillo, G. A., Israel, G. L., Belfiore, A., et al. 2020, ApJ, 895, 60
  * Sathyaprakash et al. (2019) Sathyaprakash, R., Roberts, T. P., Walton, D. J., et al. 2019, MNRAS, 488, L35
  * Swartz et al. (2011) Swartz, D. A., Soria, R., Tennant, A. F., & Yukita, M. 2011, ApJ, 741, 49
  * Titarchuk (1994) Titarchuk, L. 1994, ApJ, 434, 570
  * Tsygankov et al. (2017) Tsygankov, S. S., Doroshenko, V., Lutovinov, A. A., Mushtukov, A. A., & Poutanen, J. 2017, A&A, 605, A39
  * Vasilopoulos et al. (2021) Vasilopoulos, G., Koliopanos, F., Haberl, F., et al. 2021, ApJ, 909, 50
  * Vasilopoulos et al. (2020a) Vasilopoulos, G., Lander, S. K., Koliopanos, F., & Bailyn, C. D. 2020a, MNRAS, 491, 4949
  * Vasilopoulos et al. (2020b) Vasilopoulos, G., Ray, P. S., Gendreau, K. C., et al. 2020b, MNRAS, 494, 5350
  * Walton et al. (2011) Walton, D. J., Roberts, T. P., Mateos, S., & Heard, V. 2011, MNRAS, 416, 1844
  * Walton et al. (2015) Walton, D. J., Harrison, F. A., Bachetti, M., et al. 2015, ApJ, 799, 122
  * Walton et al. (2016) Walton, D. J., Fürst, F., Bachetti, M., et al. 2016, ApJ, 827, L13
  * Walton et al. (2018a) Walton, D. J., Fürst, F., Heida, M., et al. 2018a, ApJ, 856, 128
  * Walton et al. (2018b) Walton, D. J., Fürst, F., Harrison, F. A., et al. 2018b, MNRAS, 473, 4360
  * Wilson-Hodge et al. (2018) Wilson-Hodge, C. A., Malacaria, C., Jenke, P. A., et al. 2018, ApJ, 863, 9
  * Yoshida et al. (2010) Yoshida, T., Ebisawa, K., Matsushita, K., Tsujimoto, M., & Kawaguchi, T. 2010, ApJ, 722, 760


# SwiftAvatar: Efficient Auto-Creation of Parameterized Stylized Character
on Arbitrary Avatar Engines

Shizun Wang1, Weihong Zeng2*, Xu Wang2, Hao Yang2, Li Chen2,
Yi Yuan2, Yunzhao Zeng2, Min Zheng2 Chuang Zhang1, Ming Wu1 Corresponding
AuthorCorresponding Author

###### Abstract

††footnotetext:  *Equal Contribution††footnotetext:  $\dagger$Corresponding
Author††footnotetext: Copyright © 2023, Association for the Advancement of
Artificial Intelligence (www.aaai.org). All rights reserved.

The creation of a parameterized stylized character involves careful selection
of numerous parameters, also known as the “avatar vectors” that can be
interpreted by the avatar engine. Existing unsupervised avatar vector
estimation methods that auto-create avatars for users, however, often fail to
work because of the domain gap between realistic faces and stylized avatar
images. To this end, we propose _SwiftAvatar_ , a novel avatar auto-creation
framework that is evidently superior to previous works. SwiftAvatar introduces
dual-domain generators to create pairs of realistic faces and avatar images
using shared latent codes. The latent codes can then be bridged with the
avatar vectors as pairs, by performing GAN inversion on the avatar images
rendered from the engine using avatar vectors. Through this way, we are able
to synthesize paired data in high-quality as many as possible, consisting of
avatar vectors and their corresponding realistic faces. We also propose
semantic augmentation to improve the diversity of synthesis. Finally, a light-
weight avatar vector estimator is trained on the synthetic pairs to implement
efficient auto-creation. Our experiments demonstrate the effectiveness and
efficiency of SwiftAvatar on two different avatar engines. The superiority and
advantageous flexibility of SwiftAvatar are also verified in both subjective
and objective evaluations.

Figure 1: a) Given user’s front-facing selfie images, SwiftAvatar is able to
auto-create corresponding stylized avatars on any arbitrary avatar engine, for
example, b) the TikTok avatar engine 1, or c) the Alter avatar engine 2.

##  1 Introduction

††footnotetext:  1 https://newsroom.tiktok.com/en-us/express-yourself-through-
tiktok-avatars††footnotetext:  2 https://github.com/facemoji/alter-core

The emerging of the Metaverse concept is alongside with the wide usage of
virtual avatars, embodiment of Metaverse users in various styles, which are
popular in modern digital lives such as socialization, e-shopping and gaming.
Even though many avatar platforms, such as Zepeto, BitMoji and ReadyPlayerMe,
have enabled users to create their own stylized avatars by specifying avatar
vectors that can be interpreted by their avatar engines, the manual creation
process is tiresome and time consuming, especially when the engine provides
with a large set of options making up a long avatar vector. Besides, the
avatar vectors to specify usually consists of parameters in both continuous
forms (to control facial shape, eye spacing, etc.) and discrete forms (to
control hair styles, wearings etc.), all requiring careful selection and
cautious adjustment to achieve a satisfactory result. Therefore, it is
valuable to study _how to automatically create a stylized avatar that best
matches the user’s input selfie_.

A straightforward way to realize avatar auto-creation is using supervised
learning: a network is trained on labeled data to predict avatar vectors based
on real face inputs. However, this requires large amount of data collection
and manual labeling, which is laborious, expensive, and is not generalizable
across engines. Because the definition of avatar vectors and assets vary from
engine to engine, data labeled for one engine can not be used to train on
other engines.

Several unsupervised learning methods have been proposed to address avatar
auto-creation without using any labeled data, including Tied Output Synthesis
(TOS) (Wolf, Taigman, and Polyak 2017) and the Face-to-Parameter (F2P) series
(Shi et al. 2019, 2020). The main idea of these works can be abstracted by
Fig. 2-a. In order to achieve an avatar vector that renders avatar image as
similar to the input face as possible, these methods impose constraints on the
image-level. They suffer from issues as follows: 1) The image-level similarity
constraints they establish are designed for realistic avatar images, not
applicable to stylized avatar images that have domain gap with real face
images. 2) Leveraging the image-level supervision requires that the avatar
rendering process is differentiable. So they usually introduce an imitator
network that imitates the behavior of the non-differentiable avatar engine.
However, the un-avoidable deviation of imitators from the original avatar
engine, as illustrated by Fig. 2-a, degrades the accuracy of the similarity
measure. 3) Some approaches like F2P (Shi et al. 2019) need iterative
optimization to guarantee the quality of estimated avatar vectors, which is
time-consuming in inference.

To address these issues, we propose _SwiftAvatar_ , a novel avatar auto-
creation framework shown by Figure. 2-b). Unlike previous works that use
similarity constraints to find the avatar vector whose rendered image best
matches a given face, the core idea of our framework is cross-domain data
synthesis. SwiftAvatar is able to synthesize pairs of avatar vectors and
corresponding realistic faces in high fidelity as many as possible. They are
used to train a light-weight estimator that directly predicts avatar vectors
from input selfie. Specifically, the SwiftAvatar framework consists of three
components: _dual-domain generators_ , a pipeline for _cross-domain paired
data production_ , and an _avatar estimator_.

The dual-domain generators comprise a realistic generator and an avatar
generator, both adopting an architecture from SemanticStyleGAN (Shi et al.
2022). The realistic generator is pretrained on real faces, but the avatar
generator is transfer-learned on engine-rendered avatar images with color
consistency constraints. So that given a same latent code, the two generators
could generate a realistic face image and an avatar image that naturally look
similar. Data synthesis relies on the dual-domain generators. The production
of synthetic data starts from randomly sampled avatar vectors. They are sent
to the avatar engine to render avatar images, which are then inverted into
latent codes through the avatar generator. Finally, the latent codes are fed
into the realistic generator to get realistic face images corresponding to the
sampled avatar vectors. Moreover, we introduce semantic augmentation to expand
the diversity of produced data by adding local perturbation to latent code.
With the synthetic data, we can then train an avatar estimator on them. The
avatar estimator is light-weight and efficient at inference. Given user selfie
image, it accomplish the auto-creation by directly predicting an avatar vector
that matches the input.

Our experiments involve both objective and subjective evaluations, comparing
SwiftAvatar with previous methods in different aspects. SwiftAvatar achieves
advantageous results over all existing methods in creating avatars fidelity to
the input images. Moreover, experiments on two diverse avatar engines verify
the strong generality of SwiftAvatar. Qualitative results illustrated in Fig.
1, show that SwiftAvatar can generate reasonable avatars for both engines
given input faces. In summary, our contributions are as follows:

  * •

A novel framework, SwiftAvatar, is proposed that can automatically create a
stylized avatar given a user selfie image. It can be swiftly applied to any
arbitrary avatar engines without extra assumptions (e.g. differentiable
rendering capability) on that engine.

  * •

SwiftAvatar presents a novel pipeline that produces paired data across
domains. It involves dual-domain generators to address the domain gap between
realistic faces and stylized avatars. A novel semantic augmentation is also
devised to improve the diversity of data synthesis.

  * •

Experiments show that SwiftAvatar outperforms previous methods in terms of
both quality and efficiency. Results on two different avatar engines also
verify the strong generalizability of SwiftAvatar.

Figure 2: Conceptual comparison of previous and our methods. a) Previous
methods try to solve unknown avatar vectors under image-level constraints,
which requires to train a differentiable engine imitator to participate the
optimization. b) While our method can use origin avatar engine and cross-
domain generation to synthesize paired data (image-label), and then a model is
trained based on the synthetic dataset.

##  2 Related Work

#### 3D Face Reconstruction

Many progresses on 3D face reconstruction (Tuan Tran et al. 2017; Dou, Shah,
and Kakadiaris 2017) cannot be achieved without 3D morphable models (3DMM)
(Blanz and Vetter 1999) and its variants like BFM (Gerig et al. 2018) and
FLAME (Li et al. 2017), where the geometry and texture of a photorealistic 3D
face are parameterized as a vector through linear transform. Recent researches
also explore representing 3D faces in other formats like dense landmarks (Feng
et al. 2018b) or position maps (Feng et al. 2018a). However, stylized avatar
auto-creation is not 3D face reconstruction. They differ in two folds: 1) 3D
face reconstruction methods aim to recover realistic 3D faces, not stylized
characters; 2) Most avatar engines utilize avatar vectors to render stylized
3D characters according to manually designed assets. Estimating avatar vectors
is much more difficult than reconstructing 3DMM coefficients, since the
graphics rendering of avatar engines are usually black-box and non-
differentiable.

#### GAN and GAN Inversion

The rapid growth of generative networks, such as GANs (Goodfellow et al. 2014;
Karras et al. 2017) and VAEs (Kingma and Welling 2013; Razavi, Van den Oord,
and Vinyals 2019), inspires methods to use _latent codes_ that can implicitly
parameterize face images through a pre-trained generator. Among the pre-
trained generators, the most popular one is the generator of the StyleGAN
(Karras, Laine, and Aila 2019; Karras et al. 2020), known for its impressive
capability in generating high quality face images. In order to leverage the
pre-trained StyleGAN generator in editing existing face images, GAN inversion
(Zhu et al. 2020) is required to compute the latent code whose StyleGAN output
best matches a given face image. Existing GAN inversion methods estimate
latent codes using either iterative optimization (Abdal, Qin, and Wonka 2020)
or training a feed-forward encoder (Richardson et al. 2021) using VGG
perceptual loss (Johnson, Alahi, and Fei-Fei 2016), LPIPS loss (Zhang et al.
2018) or face identity loss (Deng et al. 2019). Though GAN inversion can not
be directly applied to estimate the avatar vectors due to the hand-crafted
nature of avatar engines, SwiftAvatar leverages GAN inversion in its data
production pipleline to build the correspondence between avatar vectors and
latent codes.

#### Portrait Stylization

Recent literatures on portrait stylization also benefit a lot from pre-trained
StyleGANs. On one hand, finetuning the pre-trained StyleGAN generator provides
with an agile approach for synthesizing faces of new styles (Song et al. 2021;
Back 2021). On the other hand, freezing low-level layers of the generator
during finetuning helps preserving structural information (e.g. face
attributes) between face images generated from the original generator
(realistic domain) and the finetuned generator (stylized domain), when they
are using similar latent codes (Back 2021; Huang, Liao, and Kwong 2021).
Though this cross-domain latent code sharing strategy cannot be directly
applied on avatar engines, we discover that it has a strong potential in
cross-domain data synthesis. In specific, we design the data production
pipeline of SwiftAvatar to leverage a pre-trained SemanticStyleGAN (Shi et al.
2022). It adopts a compositional generator architecture, disentangling the
latent space into different semantic areas separately, and could provide more
precise local controls on synthesized face images.

#### Stylized Avatar Auto-Creation

Auto-creating avatars for individual users has become an important capability
to the entertainment industry.To ensure quality, commercial solution usually
involve a large amount of manual annotations, something this paper seeks to
avoid. Among the published approaches that avoid manual annotations, some (Shi
et al. 2019, 2020) are designed for realistic avatars only, which share the
same domain with the real face images, it is easy to define whether they match
or not. Therefore, these methods can utilize straightforward in-domain
supervisions to improve fidelity of creation, such as the image-level L1/L2
losses, face parsing loss, or the face identity loss (Deng et al. 2019).
Creating stylized avatars from real faces, on the contrary, is more difficult
than creating realistic avatars. There are only few works toward this
direction. Tied Output Synthesis (TOS) (Wolf, Taigman, and Polyak 2017)
devises a encoder-decoder structure, that eliminates the domain gap by sharing
one encoder across two domains. The contemporaneous work AgileAvatar (Sang et
al. 2022) formulates a cascaded framework which progressively bridges the
domain gap. To address the non-differentiable issue of avatar engines, all
these methods need to train an imitator network to imitate the engines’
rendering procedure. While the quality of the images generated by the imitator
is relatively poor. Besides, they all impose explicit constraints on images
from different domains, resulting in suboptimal matches because of the domain
gap. By contrast, our SwiftAvatar framework employs original avatar engine to
generate high-quality images and utilizes dual-domain generators to overcome
the domain gap.

Figure 3: Overview of our method. a) Dual-domain generators consist of a fixed
realistic generator and a transfer-learned avatar generator. They can generate
corresponding realistic and avatar images when given the same latent codes. b)
Paired data production starts from randomly sampled avatar vectors, which are
sent to the avatar engine to render avatar images. These images are then
inverted into latent codes through GAN inversion, and the latent codes are fed
into the realistic generator to get realistic face images corresponding to the
sampled avatar vectors. c) Supervised by the paired data, the avatar estimator
can predict avatar vectors when given an real face input, then the avatar can
be rendered by avatar engine.

##  3 Methodology

In this section, we present our unsupervised framework for stylized avatar’s
auto-creation, as shown in Fig.3. It aims at estimating avatar vector $p$ for
input real face image $x$, then the estimated $p$ can be used to render the
corresponding stylized avatar $y$ by avatar engine $E$. Our solution is split
into three parts: dual-domain generators in Sec. 3.1, paired data production
in Sec. 3.2, and avatar estimator in Sec. 3.3. Dual-domain generators address
the domain gap problem by generating realistic faces and stylized avatars with
shared latent codes. Then, paired data production focuses on how to generate
paired data consisting of avatar vectors and realistic faces. Finally, avatar
estimator estimates desired avatar vectors similar to input real face images.

###  3.1 Dual-domain Generators

The dual-domain generators consist of a realistic generator
$G_{\textit{real}}$ and an avatar generator $G_{\textit{avatar}}$ to achieve
cross-domain image generation. Given the same latent code, they can
simultaneously synthesize paired images of both realistic face and stylized
avatar while preserving the same attributes (e.g. skin color, hair style). To
impose an extra facial consistency constraint between the realistic and the
avatar domains, we adopt SemanticStyleGAN(Shi et al. 2022) as the architecture
of two generators owing to its extra semantic segmentation output.

#### Cross-domain Generation

Pretrained SemanticStyleGAN on CelebAMask-HQ (Lee et al. 2020) is directly
used as realistic generator $G_{\textit{real}}$, and also used to initialize
the weight of $G_{\textit{avatar}}$. We perform transfer learning on
$G_{\textit{avatar}}$: using only limited number of avatar images
$\mathcal{Y}$ to finetune avatar generator $G_{\textit{avatar}}$.
$\mathcal{Y}$ are rendered from avatar engine using randomly sampled avatar
vectors. The finetune procedure follows the settings in SemanticStyleGAN,
which using the loss of StyleGAN2 (Karras et al. 2020):

$\mathcal{L}_{StyleGAN2}=\mathcal{L}_{adv}+\lambda_{R1}\mathcal{L}_{R1}+\lambda_{path}\mathcal{L}_{path}$
(1)

where $\lambda_{R1}$, $\lambda_{path}$ are the constant weights of R1
regularization (Mescheder, Geiger, and Nowozin 2018) and path length
regularization (Karras et al. 2020) separately. Adversarial loss
$\mathcal{L}_{adv}$ adopts non-saturating logistic loss (Goodfellow et al.
2014) and it forces $G_{\textit{avatar}}$ to generate images similar to avatar
image dataset $\mathcal{Y}$. R1 regularization $\mathcal{L}_{R1}$ is employed
to improve the training stability and reduce the number of artifacts. And path
length regularization $\mathcal{L}_{path}$ leads to more reliable and
consistently behaving models.

#### Facial Consistency

Although directly fine-tuning SemanticStyleGAN can generate structurally
similar paired data of avatar image and realistic face, the colors of each
region are not well matched. Since SemanticStyleGAN learns the joint modeling
of image and semantic segmentation. Such design can simultaneously synthesize
face images and their semantic segmentation results. We utilize the semantic
segmentation output and introduce a color matching loss for cross-domain
facial color consistency. Specifically, we extract specified pixels from same
semantic areas in generated paired images, and match the mean color of them.
$m^{s}(I)$ is the mean color of the region $s$ in the image $I$, and we mainly
consider matching the color in hair and skin areas: $s\in\\{hair,skin\\}$. The
color matching loss is:

$\mathcal{L}_{color}=\sum_{s}\left\|m^{s}(G_{\textit{real}}(z))-m^{s}(G_{\textit{avatar}}(z))\right\|^{2}$
(2)

Overall, the final finetuning loss $\mathcal{L}_{total}$ for
$G_{\textit{avatar}}$ is formulated as:

$\mathcal{L}_{total}=\mathcal{L}_{StyleGAN2}+\lambda_{color}\mathcal{L}_{color}$
(3)

An example is shown in Figure. 4, when given a shared latent code, the dual-
domain generators can generate a pair of images: a realistic face and a
stylized avatar. They are similar in facial structure and color composition,
yet belong to two different domains.

###  3.2 Paired Data Production

Paired data production pipeline focuses on synthesizing paired avatar vectors
and realistic face images, as is illustrated in Figure. 3. We sample a number
of random avatar vectors as labels $\mathcal{P}$, which are used by the
graphics engine to generate corresponding avatar images $\mathcal{Y}$. Then,
for every avatar image $y$, we pass it through GAN inversion to get its latent
code $w$. We adopt optimization-based GAN inversion for its better
performance:

$w^{*}=\mathop{\arg\min}\limits_{w}\mathcal{L}_{invert}(y,w)$ (4)

where $w^{*}$ represents our target latent code. For faster convergence,
latent codes are initialized with mean value $w_{mean}$, and optimized by
gradient descent. We use LPIPS loss (Zhang et al. 2018) to measure perceptual
similarity, and mean squared error (MSE) loss between original images and
reconstructed avatar images to measure reconstruction similarity. Besides, MSE
loss between $w_{mean}$ and $w$ is set as latent regularization to ensure
better generation results. The loss function is formulated as:

$\displaystyle\mathcal{L}_{invert}=\ $ $\displaystyle\lambda_{p}\
\textit{LPIPS}(y,G_{\textit{avatar}}(w))+$
$\displaystyle\lambda_{i}\left\|y-G_{\textit{avatar}}(w)\right\|^{2}+\lambda_{l}\left\|w-w_{mean}\right\|^{2}$
(5)

where $\lambda_{p}$, $\lambda_{i}$, $\lambda_{l}$ are constants to balance
three loss terms. After obtaining desired latent code $w$, it is fed into the
realistic generator $G_{\textit{real}}$ to generate a realistic face $x$ which
is similar to original avatar image $y$ in identification:

$x=G_{\textit{real}}(w)$ (6)

In this way, we can generate paired data $(p,x)$ used for later avatar vector
estimation training process. An example illustration of paired data is shown
in Figure. 5.

Figure 4: Generated results of our dual-domain generators. When providing
randomly sampled latent codes, we can obtain a) and b), realistic faces and
their semantic segmentation results. Simultaneously, we can obtain c) and d),
stylized avatar images and their semantic segmentation results.

#### Semantic Augmentation

The sampled avatar vectors, as well as their rendered avatar faces suffer from
the lack of diversity, due to the limited amount of assets available for
avatar engines. Take the “hair style” attribute for example, in most stylized
avatar engines, different hair styles are determined by selecting from
different hair meshes all predefined in the engine. The limited number of
predefined mesh assets restricts the capability of avatar engines to match a
real-world user face, whose hair styles would be countless. To enrich the
diversity of generated realistic faces in data production, we take advantage
of the compositional generation ability of our dual-domain generators which
are implemented as SemanticStyleGANs, and design the _semantic augmentation_.
In SemanticStyleGAN, each semantic part is modulated individually with
corresponding local latent codes. Such property enables us to manipulate only
specific regions of the synthetic images while keeping other regions
unchanged. In implementation, we add random noises to part of latent code
corresponding to these ambiguous avatar vectors (e.g. hair type). Semantic
augmentation can be described as:

$w_{k}\to(1-\lambda_{aug})w_{k}+\lambda_{aug}n$ (7)

where $\lambda_{aug}$ is a hyper-parameter to adjust semantic augmentation
tensity, $w_{k}$ is the local part of latent code to be augmented, and $n$
represents for random noise. A semantic augmentation example is shown in
Figure. 5, where hair and background region are changed.

Figure 5: Examples of produced paired data. a) Engine rendered avatar images.
b) Reconstructed avatar images. c) Generated realistic images. d) Semantic
augmented images.

###  3.3 Avatar Estimator

Once the aforementioned synthetic paired dataset is produced, we can train an
avatar estimator to predict avatar vectors, which contain continuous and
discrete parameters. We choose ResNet-18 (He et al. 2016) pretrained on
ImageNet (Deng et al. 2009) as our backbone. We remove its last fully
connected layers and add multiple separate MLP heads for different parameter
estimation. All continuous parameters form a target vector to be predicted in
one head, supervised by L1 loss. Every discrete parameter estimation is
carried out with a standalone head. Because both generation and semantic
augmentation would inevitably introduce noises to discrete labels, we choose
symmetric cross-entropy (SCE) loss (Wang et al. 2019), which has been proven
robust to noises, for the optimization of discrete tasks. The total loss of
avatar estimator is:

$\mathcal{L}_{estimator}=\lambda_{d}\
\textit{SCE}(\widehat{p_{d}^{i}},p_{d}^{i})+\lambda_{c}\left|\widehat{p_{c}}-p_{c}\right|$
(8)

where $\lambda_{d}$ and $\lambda_{c}$ are hyper-parameters to balance two loss
terms. $\widehat{p_{d}^{i}}$, $\widehat{p_{c}}$ are the prediction results of
$i$-th discrete head and continuous head respectively. And ${p_{d}^{i}}$,
${p_{c}}$ are their corresponding ground-truth.

Figure 6: Visual comparison with other methods on TikTok engine and Alter
engine. a) Given input images, b) our method creates stylized avatars that
look similar to input images. c) Baseline method without considering the
domain gap problem in stylized avatars. d) F2P (Shi et al. 2019), an image-
level constraint method intending to create semi-realistic avatars, is not
suitable to create stylized avatars. e) F2P v2 (Shi et al. 2020), a fast and
robust version of F2P. f) Manually created avatars by professional designers,
which can be regarded as ground truth.

##  4 Experiments

#### Experimental Data

To verify the effectiveness of our method, we conduct experiments on two
stylized avatar engines: the TikTok engine and the Alter engine. The TikTok
engine contains resourceful discrete and continuous avatar parameters, so we
generate 50000 avatar vectors and corresponding rendered images. The Alter
engine is an open-source avatar engine and only contains discrete avatar
parameters which has fewer assets than the TikTok engine, so we just generate
10000 avatar vectors and corresponding rendered images. The detailed
information of both avatar engines can be found in supplementary materials.
These avatar images and avatar vectors are used for finetuning avatar
generator and producing paired data. For evaluation, we choose 116 images from
FFHQ dataset (Karras, Laine, and Aila 2019), which consists of diverse kinds
of face shapes, hairstyles, etc. We invite designers to manually generate
avatar vectors for these 116 images as ground truth.

#### Implementation Details

We implement our methods using PyTorch 1.10 library and perform all
experiments on NVIDIA V100 GPUs. When finetuning the avatar generator
$G_{\textit{avatar}}$, we use the same optimizer settings as in
SemanticStyleGAN. Batch size is set to 16, style mixing probability (Karras,
Laine, and Aila 2019) is set to 0.3. $\lambda_{R1}$, $\lambda_{path}$ are set
to 10 and 0.5 separately. Lazy regularization (Karras et al. 2020) is applied
every 16 mini-batches for discriminator (R1 regularization) and every 4 mini-
batches for generator (path length regularization). All the images used for
generators are aligned and resized to resolution $512\times 512$. The
optimization-based GAN inversion approach employs Adam (Kingma and Ba 2014)
optimizer in the paired data production stage, and the learning rate initially
follows cosine annealing with 0.1. We optimize 200 steps for all latent codes,
and $\lambda_{i}$, $\lambda_{p}$, $\lambda_{l}$ are set to 0.1, 1 and 1,
respectively. The mean latent code $w_{mean}$ is the average among $10^{4}$
randomly generated latent codes from avatar generator $G_{\textit{avatar}}$ in
$\mathcal{W}$ space, and serve as the initialization. Notably, directly
optimizing the latent code could be problematic since some avatar assets are
transparent, e.g. glasses. Thus, we use a modified version $\tilde{w}_{mean}$
for latent code initialization (See supplementary materials for details). For
semantic augmentation, we generate 10 augmented images for each latent code,
using randomly generated noise in $\mathcal{W}$ space.we set $\lambda_{aug}$
to 1 for the background to improve the model robustness of background
variance, and also set $\lambda_{aug}$ to 0.3, 0.06 for the hair part and
glasses part to expand data diversity. In the avatar estimator training stage,
the input images of avatar estimator are resized to $224\times 224$. We use
the Adam optimizer with batch size 256 to train 100 epochs. The learning rate
is set to $1e-3$, and decayed by half per 30 epochs. For the experiments on
TikTok engine, there are 1 continuous head and 8 discrete heads inside the
avatar estimator, so we set $\lambda_{c}$ and $\lambda_{d}$ to 1 and 10
separately. For the experiments on Alter engine, the avatar estimator contains
6 discrete heads and the training loss only contains discrete loss.

###  4.1 Comparison with Other methods

We compare our method with other methods, including Baseline, F2P (Shi et al.
2019), F2P v2 (Shi et al. 2020) on both their auto-creation similarity and
inference speed. Baseline method is setup as ignoring the domain gap problem,
where the avatar estimator is trained on $(y,p)$ paired data, that is, trained
on rendered avatar images instead of real face images. Figure. 6 shows a
comparison of rendered stylized avatars among different methods corresponding
to their predicted avatar vectors. As can be seen, our method can better
address the domain gap problem and create stylized avatars with high
similarity with input real faces and approximate the quality of manual method
from designers. For more results, please refer to supplementary materials.

Method | TikTok engine $\downarrow$ | Alter engine $\downarrow$ | Speed $\uparrow$
---|---|---|---
Baseline | 0.5033 | 0.3962 |  $\sim 10^{2}$ Hz
F2P | 0.3562 | 0.3040 |  $\sim 1$ Hz
F2P v2 | 0.4466 | 0.2825 |  $\sim 10^{2}$ Hz
Ours | 0.3110 | 0.2405 |  $\mathbf{\sim 10^{2}}$ Hz
Table 1: Quantitative Evaluation. We compare our method with three other
methods on TikTok engine and Alter engine. Lower distance represents results
are more similar to manual-creation. We test speed on NVIDIA V100.

###  4.2 Quantitative Evaluation

Although the avatar auto-creation has no standard answer, we generally regard
the avatar manually created by professional designers as the ground truth, and
quantitatively evaluate at image level. We calculate the perceptual distance
(LPIPS) (Huang, Liao, and Kwong 2021) between auto-created avatar and manual-
created avatar to simulate human observation. The lower distance indicates the
avatar is better matched to input image. The results are presented in Table.
1, from which we can see our method significantly outperforms others. Since
our method only needs one forward propagation, our inference speed is also
competitive, which can be applied in real-time applications.

###  4.3 Human Rating

We invited 50 volunteers to subjectively evaluate all algorithms of 50 images
from evaluation dataset on both graphics engine. Given the input real face
images and stylized avatars generated from our and other methods, the
volunteers were requested to pick the best matching avatar for each input real
face. The methods were presented in random order, and the volunteers were
given an unlimited time to choose. The human rating results are shown in
Table. 2. $59.11\%$ of the answers on TikTok engine and $63.68\%$ on the Alter
engine selected our method as the best matching results, showing the
superiority of our method compared with others.

Method | Ours | baseline | F2P | F2P v2
---|---|---|---|---
TikTok |  59.11% | 5.66% | 30.86% | 4.37%
Alter |  63.68% | 18.94% | 8.65% | 8.73%
Table 2: Human subjective rating on two engines. Our method earns the most choices when asked to choose the avatar which matches the human image. Method | TikTok engine $\downarrow$ | Alter engine $\downarrow$
---|---|---
baseline | 0.5033 | 0.3962
\+ domain adaptation | 0.3401 | 0.3123
\+ semantic aug | 0.3110 | 0.2405
Table 3: Ablation study on two engines. Lower distance indicates better
matching results to manual-creation.

###  4.4 Ablation Study

We conduct the ablation experiments in terms of domain adaptation and semantic
augmentation to verify their importance in our framework. We adopt the same
evaluation metric described in Sec. 4.2. The ablation starts from Baseline
method, where the avatar estimator is trained on rendered avatar images and
avatar vectors pair. Then on the basis of baseline method, we add domain
adaptation and semantic augmentation in turn. Table. 3 shows their
quantitative results on different engines. The domain adaptation greatly
alleviates the domain gap problem in stylized avatar auto-creation,
establishes a bridge of connecting real faces and stylized avatars. The
semantic augmentation brings noticeable improvement for our method due to its
expansion of diversity to samples.

##  5 Limitations and Future Work

There are two main limitations we observed in the experiments. First, our
method occasionally predicts wrong color influenced by environmental lighting.
This problem might be resolved by considering lighting condition into our
pipeline. Second, in paired data production, avatar vector sampling
distribution directly influences the training data quality. Simply random
sampling could produce some strange images and may cause long-tail problem
(e.g. gender, age). In the future, we will perform attribute analysis, and
introduce reasonable sampling priors to address synthetic data distribution
problem.

##  6 Conclusion

In summary, we present a novel unsupervised framework for auto-creation of
stylized avatars. We design the dual-domain generators to address the domain
gap between the real images and stylized avatars. Then following the paired
data production pipeline, high-quality paired data are produced, which is used
for training the avatar estimator. Finally stylized avatars are created by
conducting efficient avatar vector estimation. Compared with previous methods,
our method is more concise in training stage and more efficient in inference
stage. Results on quantitative evaluation and human rating demonstrate the
superiority of our method. Also, the success of applying on two different
avatar graphics engines demonstrates the generality of our method.

## References

  * Abdal, Qin, and Wonka (2020) Abdal, R.; Qin, Y.; and Wonka, P. 2020.  Image2stylegan++: How to edit the embedded images?  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 8296–8305.
  * Back (2021) Back, J. 2021.  Fine-Tuning StyleGAN2 For Cartoon Face Generation.  _arXiv preprint arXiv:2106.12445_.
  * Blanz and Vetter (1999) Blanz, V.; and Vetter, T. 1999.  A morphable model for the synthesis of 3D faces.  In _Proceedings of the 26th annual conference on Computer graphics and interactive techniques_ , 187–194.
  * Deng et al. (2019) Deng, J.; Guo, J.; Xue, N.; and Zafeiriou, S. 2019.  Arcface: Additive angular margin loss for deep face recognition.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 4690–4699.
  * Dou, Shah, and Kakadiaris (2017) Dou, P.; Shah, S. K.; and Kakadiaris, I. A. 2017.  End-to-end 3D face reconstruction with deep neural networks.  In _proceedings of the IEEE conference on computer vision and pattern recognition_ , 5908–5917.
  * Feng et al. (2018a) Feng, Y.; Wu, F.; Shao, X.; Wang, Y.; and Zhou, X. 2018a.  Joint 3d face reconstruction and dense alignment with position map regression network.  In _Proceedings of the European conference on computer vision (ECCV)_ , 534–551.
  * Feng et al. (2018b) Feng, Z.-H.; Huber, P.; Kittler, J.; Hancock, P.; Wu, X.-J.; Zhao, Q.; Koppen, P.; and Rätsch, M. 2018b.  Evaluation of dense 3D reconstruction from 2D face images in the wild.  In _2018 13th IEEE International Conference on Automatic Face & Gesture Recognition (FG 2018)_, 780–786. IEEE.
  * Gerig et al. (2018) Gerig, T.; Morel-Forster, A.; Blumer, C.; Egger, B.; Luthi, M.; Schönborn, S.; and Vetter, T. 2018.  Morphable face models-an open framework.  In _2018 13th IEEE International Conference on Automatic Face & Gesture Recognition (FG 2018)_, 75–82. IEEE.
  * Goodfellow et al. (2014) Goodfellow, I.; Pouget-Abadie, J.; Mirza, M.; Xu, B.; Warde-Farley, D.; Ozair, S.; Courville, A.; and Bengio, Y. 2014.  Generative adversarial nets.  _Advances in neural information processing systems_ , 27.
  * Huang, Liao, and Kwong (2021) Huang, J.; Liao, J.; and Kwong, S. 2021.  Unsupervised image-to-image translation via pre-trained stylegan2 network.  _IEEE Transactions on Multimedia_ , 24: 1435–1448.
  * Johnson, Alahi, and Fei-Fei (2016) Johnson, J.; Alahi, A.; and Fei-Fei, L. 2016.  Perceptual Losses for Real-Time Style Transfer and Super-Resolution.  In Leibe, B.; Matas, J.; Sebe, N.; and Welling, M., eds., _Computer Vision – ECCV 2016_ , 694–711. Cham: Springer International Publishing.  ISBN 978-3-319-46475-6.
  * Karras et al. (2017) Karras, T.; Aila, T.; Laine, S.; and Lehtinen, J. 2017.  Progressive growing of gans for improved quality, stability, and variation.  _arXiv preprint arXiv:1710.10196_.
  * Karras, Laine, and Aila (2019) Karras, T.; Laine, S.; and Aila, T. 2019.  A style-based generator architecture for generative adversarial networks.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 4401–4410.
  * Karras et al. (2020) Karras, T.; Laine, S.; Aittala, M.; Hellsten, J.; Lehtinen, J.; and Aila, T. 2020\.  Analyzing and improving the image quality of stylegan.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 8110–8119.
  * Kingma and Ba (2014) Kingma, D. P.; and Ba, J. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_.
  * Kingma and Welling (2013) Kingma, D. P.; and Welling, M. 2013.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_.
  * Lee et al. (2020) Lee, C.-H.; Liu, Z.; Wu, L.; and Luo, P. 2020.  MaskGAN: Towards Diverse and Interactive Facial Image Manipulation.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_.
  * Li et al. (2017) Li, T.; Bolkart, T.; Black, M. J.; Li, H.; and Romero, J. 2017.  Learning a model of facial shape and expression from 4D scans.  _ACM Transactions on Graphics, (Proc. SIGGRAPH Asia)_ , 36(6): 194:1–194:17.
  * Mescheder, Geiger, and Nowozin (2018) Mescheder, L.; Geiger, A.; and Nowozin, S. 2018.  Which training methods for GANs do actually converge?  In _International conference on machine learning_ , 3481–3490. PMLR.
  * Razavi, Van den Oord, and Vinyals (2019) Razavi, A.; Van den Oord, A.; and Vinyals, O. 2019.  Generating diverse high-fidelity images with vq-vae-2.  _Advances in neural information processing systems_ , 32.
  * Richardson et al. (2021) Richardson, E.; Alaluf, Y.; Patashnik, O.; Nitzan, Y.; Azar, Y.; Shapiro, S.; and Cohen-Or, D. 2021.  Encoding in style: a stylegan encoder for image-to-image translation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2287–2296.
  * Sang et al. (2022) Sang, S.; Zhi, T.; Song, G.; Liu, M.; Lai, C.; Liu, J.; Wen, X.; Davis, J.; and Luo, L. 2022.  AgileAvatar: Stylized 3D Avatar Creation via Cascaded Domain Bridging.  _arXiv preprint arXiv:2211.07818_.
  * Shi et al. (2019) Shi, T.; Yuan, Y.; Fan, C.; Zou, Z.; Shi, Z.; and Liu, Y. 2019.  Face-to-parameter translation for game character auto-creation.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 161–170.
  * Shi et al. (2020) Shi, T.; Zuo, Z.; Yuan, Y.; and Fan, C. 2020.  Fast and robust face-to-parameter translation for game character auto-creation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 1733–1740.
  * Shi et al. (2022) Shi, Y.; Yang, X.; Wan, Y.; and Shen, X. 2022.  SemanticStyleGAN: Learning Compositional Generative Priors for Controllable Image Synthesis and Editing.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 11254–11264.
  * Song et al. (2021) Song, G.; Luo, L.; Liu, J.; Ma, W.-C.; Lai, C.; Zheng, C.; and Cham, T.-J. 2021\.  AgileGAN: stylizing portraits by inversion-consistent transfer learning.  _ACM Transactions on Graphics (TOG)_ , 40(4): 1–13.
  * Tuan Tran et al. (2017) Tuan Tran, A.; Hassner, T.; Masi, I.; and Medioni, G. 2017.  Regressing robust and discriminative 3D morphable models with a very deep neural network.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 5163–5172.
  * Wang et al. (2019) Wang, Y.; Ma, X.; Chen, Z.; Luo, Y.; Yi, J.; and Bailey, J. 2019.  Symmetric cross entropy for robust learning with noisy labels.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 322–330.
  * Wolf, Taigman, and Polyak (2017) Wolf, L.; Taigman, Y.; and Polyak, A. 2017.  Unsupervised creation of parameterized avatars.  In _Proceedings of the IEEE International Conference on Computer Vision_ , 1530–1538.
  * Zhang et al. (2018) Zhang, R.; Isola, P.; Efros, A. A.; Shechtman, E.; and Wang, O. 2018.  The unreasonable effectiveness of deep features as a perceptual metric.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 586–595.
  * Zhu et al. (2020) Zhu, J.; Shen, Y.; Zhao, D.; and Zhou, B. 2020.  In-domain gan inversion for real image editing.  In _European conference on computer vision_ , 592–608. Springer.

## Appendix

### Implementation Details

In the process of GAN-inversion, we optimize the latent code with
initialization $w_{mean}$. Notably, As glasses is an optional discrete avatar
vector, normally using $w_{mean}$ may cause the absence of glasses when it
should appear due to their transparent design in SemanticStyleGAN. Since the
latent code of SemanticStyleGAN is disentangle to different local part, we
select the latent codes which contains the glasses and average them to
$w_{mean}^{glasses}$, then extract the glasses part and replace it into
$w_{mean}$ to form a new mean latent code $\tilde{w}_{mean}$, it can generate
glasses correctly after optimization and we use it in practical.

Components | Vectors | Numbers | Value Type
---|---|---|---
head | head type | 6 | continuous
head | skin tone | 39 | discrete
head | head width | 1 | continuous
head | head length | 1 | continuous
eye | eye type | 4 | continuous
eye | eye rotation | 1 | continuous
eye | eye spacing | 1 | continuous
eye | eye size | 1 | continuous
mouth | mouth type | 9 | discrete
mouth | mouth width | 1 | continuous
mouth | mouth volume | 1 | continuous
mouth | mouth position | 1 | continuous
nose | nose type | 5 | discrete
nose | nose width | 1 | continuous
nose | nose height | 1 | continuous
nose | nose position | 1 | continuous
hair | hair type | 36 | discrete
hair | hair color | 9 | discrete
brow | type | 8 | discrete
facial hair∗ | type | 4 | discrete
glasses∗ | type | 1 | discrete
Table 4: A detailed description of avatar vector in TikTok engine. (*
indicates optional)

### Avatar Engine Information

Here we provide detailed information of avatar engine we used in the
experiments, which are shown in Table. 4 and Table. 5. ”Component” represents
avatar parts which avatar vectors belong to. ”Vectors” represents the
parameter users can adjust (continuous) or select (discrete). ”Numbers” of
discrete vectors represents the assets numbers, ”Number” of continuous vectors
is 1. ”Value Type” indicates whether continuous or discrete vectors.

The information of TikTok avatar engine is listed in Table. 4. It contains 8
discrete avatar vectors and 13 continuous vectors. And Alter avatar engine
only contains 6 discrete avatar vectors, which is listed in Table. 5. The
experiments on two diverse avatar engines demonstrates the generality of our
method.

Components | Controllers | numbers | value type
---|---|---|---
head | head type | 6 | discrete
head | skin tone | 8 | discrete
hair | hair type | 22 | discrete
hair | hair color | 7 | discrete
brow | type | 5 | discrete
glasses∗ | type | 1 | discrete
Table 5: A detailed description of avatar vector in Alter engine. (* indicates
optional)

### More Visual Results

Figure. 7 and Figure. 8 show more visual comparison on TikTok engine and Alter
engine.

Figure 7: More results on TikTok engine: Input, Ours, Baseline, F2P, F2P v2,
Manual Figure 8: More results on Alter engine: Input, Ours, Baseline, F2P, F2P
v2, Manual