Unveiling the origin of the optical/UV emission from the Galactic ULX Swift J0243.6+6124 during its 2017-2018 giant outburst

From late September 2017 to February 2018, the Be X-ray binary (BeXB) Swift J0243.6+6124 underwent an unprecedently bright giant outburst. The reported X-ray luminosities were so high that the system was classified as an Ultraluminous X-ray source (ULX). It was also the first BeXB pulsar showing radio jet emission. The source was not only bright in X-rays and radio, but also in optical and UV wavelenghts. In this paper we aim to understand the origin of the observed optical/UV fluxes simultaneous to the X-ray emission. We studied the optical/UV light curves in comparison with the X-ray fluxes along the outburst, considering the main mechanisms that can explain the optical/UV emission in X-ray binaries. Due to the tight correlation observed between the optical/UV and X-ray light curves, reprocessing of X-rays seems to be the most plausible explanation. We calculated the timescales of the light curves decays and studied the correlation indexes between the optical and X-ray emission. Finally, we built a physical model considering X-ray heating of the surface of the donor star, irradiation of the accretion disk, and emission from a viscously heated accretion disk, in order to reproduce the observed optical/UV SEDs along the outburst. We considered in our model that the Be circumstellar disk was co-planar to the orbit, and then we neglected its irradiation in the current model. As an input of the model, we used as incident X-ray luminosities those calculated from the bolometric X-ray fluxes obtained from the spectral fit of the Swift/XRT and BAT observations. We conclude that reprocessing of X-rays as X-ray heating of the Be star surface and irradiation of the accretion disk are the two main mechanisms that can reproduce the observed optical/UV emission during the 2017-2018 giant outburst of Swift J0243.6+6124.


Introduction
Be/X-ray binaries (BeXBs) are systems formed by a neutron star (NS) and a Be star.The Be star expels matter from the pho-(L X ≈ 10 36 -10 37 erg s −1 ); and type II or giant outbursts, with L X > 10 37 erg s −1 .Type I or normal outbursts take place close to the periastron passage of the NS, when the NS can cross the circumstellar disk and accrete material from it emitting in X-rays.Type II outbursts can occur at any orbital phase, with a duration usually longer than the orbital period, and their origin remains unclear, although two main scenarios have been proposed to explain their ocurrence.In one hand, type II outbursts could take place due to the interaction of the NS with a tidally warped and precessing circumstellar disk (Moritani et al. 2011(Moritani et al. , 2013;;Okazaki et al. 2013).On the other hand, for highly misaligned circumstellar disks, giant outbursts could take place when the Be disk undergoes strong eccentricity growth due to Kozai-Lidov oscillations and the material overflows the Roche Lobe and is captured by the NS (Martin et al. 2014;Laplace et al. 2017;Martin & Franchini 2019).
Swift J0243.6+6124 is a unique system, since it is the first and only Ultraluminous X-ray Pulsar (ULXP) in our Galaxy (Doroshenko et al. 2018).It was discovered by the Burst Alert Telescope (BAT) on board the Neil Gehrels Swift observatory on 3 October 2017 at the beginning of a giant X-ray outburst (Kennea et al. 2017).A spin period of around 9.86 s was detected by Swift/XRT (Kennea et al. 2017), Fermi/GBM (Jenke & Wilson-Hodge 2017), and NuSTAR (Bahramian et al. 2017).The Fermi/GBM observations allowed accurate determination of the orbital period P orb =27.587 days and eccentricity e=0.098 (Wilson-Hodge et al. 2018;Zhang et al. 2019).
A few days after the discovery of the system, Kouroubatzakis et al. ( 2017) obtained an optical spectrum of the source and claimed it was a new Be/X-ray binary (BeXB).A detailed analysis of the spectroscopic and photometric observations of the optical counterpart of Swift J0243.6+6124 was performed by Reig et al. (2020) (hereafter Paper I), who classified the optical companion of this system as an O9.5Ve star.In this work, long-term optical variability was observed in the available photometric observations, which was interpreted as variations in the Be circumstellar disk.
Different estimations of the distance to the source have been obtained.Bikmaev et al. (2017) proposed a distance of ∼ 2 kpc based on optical photometric observations.Doroshenko et al. (2018) and Zhang et al. (2019) proposed a lower limit of ∼ 5 kpc based on the spin-up rate and accretion-torque models.An estimation of 4.8±0.5 kpc was obtained from optical photometric observations in Paper I. It should be stressed that many papers focused on the X-ray analysis of the 2017 giant outburst of Swift J0243.6+6124 have considered in their analyses the distance estimated from the Gaia Data Release 2 at face value, which is d=6.8 kpc (Bailer-Jones et al. 2018).However, recent estimations from the Gaia Early Data Release 3 points to a distance value of d=5.2±0.3 kpc (Bailer-Jones et al. 2021), in agreement with the estimation from the optical data in Paper I. We will use this value in our calculations in this work.This value of the distance implies X-ray luminosities at the peak of the outburst close to L X ∼ 10 39 erg/s, which still indicates that the Eddington limit for the neutron star was exceeded during the outburst, and the system can still be classified as an Ultraluminous X-ray Source (ULX).
Several studies have been focused on the analysis of the X-ray spectrum, which has been characterized by a cut-off power-law continuum (Zhang et al. 2019) with an evolving iron line emission (Jaisawal et al. 2019).Spectral and pulse profile analyses have been performed by Wilson-Hodge et al. (2018); Doroshenko et al. (2020) and Liu et al. (2022a).There have been also many studies focused on studying the configuration of the magnetic field of Swift J0243.6+6124(Doroshenko et al. 2018;Wilson-Hodge et al. 2018;Tsygankov et al. 2018;Doroshenko et al. 2020;Bykov et al. 2022;Serim et al. 2023), and a cyclotron line from 120 to 146 keV was proposed by Kong et al. (2022).One or more blackbody components associated with the hot spot around the polar region and thermal emission from the accretion column and from the photosphere of a possible ultrafast outflow have been proposed to fit the X-ray spectrum at lower energies (Tao et al. 2019;van den Eijnden et al. 2019b).
An interesting and surprising feature of the source was the detection of radio emission during the X-ray outburst, which has been interpreted as coming from a jet (van den Eijnden et al. 2018).Radio emission was not expected for any BeXB, since it was believed that strong magnetic fields should inhibit the jet formation (Fender & Hendry 2000;Massi & Kaufman Bernadó 2008).Anyway, the radio emission detected in Swift J0243.6+6124 was two orders of magnitude fainter than for other X-ray pulsars in low-mass X-ray binaries (LMXBs) at similar X-ray luminosities, and this difference could be due to the strong magnetic field.During a later X-ray re-brightening of the source at significantly lower X-ray luminosity, the radio emission re-brightened to similar levels as the main outburst peak (van den Eijnden et al. 2019a).
In Paper I, we identified a very bright optical outburst in the V-Johnson band optical light curve of Swift J0243.6+6124,taking place simultaneously with the X-ray emission.In this paper, we explore the possible mechanisms that can explain the optical and UV emission observed during the giant outburst of the system and the correlation with the X-ray emission.Some of these mechanisms involve variations in the companion star and/or the circumstellar disk of the Be star and can produce optical variability of tenths of magnitudes.They comprise different scenarios such as ejections of matter from the Be star (Hubert & Floquet 1998;Mennickent et al. 2002), precession and/or warping or changes in the inclination of the circumstellar disk (Okazaki et al. 2013;Martin et al. 2014;Martin & Franchini 2019), and formation and dissipation of the circumstellar disk (see for example Reig et al. 2023).
On the other hand, there are other mechanisms directly related to the X-ray emission.The correlation indexes between the X-ray and optical/IR fluxes have been traditionally used in the past to study the origin of the optical emission in X-ray binaries.van Paradijs & McClintock (1994) affirmed that if the optical emission comes from reprocessing of X-rays in the accretion disk of LMXBs, then L V ∝ T 2 ∝ L 0.5 X a, where T is the effective temperature of the irradiated accretion disk and a is the orbital separation.Statistical studies have identified the emission from irradiated disks as the main mechanism contributing to the optical and IR light in many LMXBs (Russell et al. 2006).For neutron star LMXBs, the origin of the optical/IR emission has been found to vary with X-ray luminosity and source type (atolls, millisecond X-ray pulsars, and Z sources).Russell et al. (2007) concluded that X-ray reprocessing can explain the optical/IR emission in Z sources at all luminosities.For atoll and millisecond pulsars this effect would dominate at low luminosities (L X ∼ 10 36 erg/s and L X ∼ 10 37 erg/s respectively), while above L X ⪆ 10 37 erg/s, optically thin synchrotron emission from the jet would dominate the NIR light (Russell et al. 2007).King & Ritter (1998) predicted that when the optical emission comes mainly from an irradiated disk, the exponential decay timescales of the optical and IR light curves should be roughly ∼2-4 times the e-folding time of the X-ray light curve.This relationship has been observed in several X-ray novae, with typical X-ray light curve decay timescales of 30 days (Chen et al. 1997), Article number, page 2 of 16 and more recently also for the black hole candidate X-ray transient XTE J181-330 (Rykoff et al. 2007).
Statistical studies are useful to study the general behaviour of X-ray binaries and to find patterns for different types and subtypes, or when the sources are at different spectral states or luminosities.However, simple flux-flux correlations can be too simplistic to understand in detail the effects of reprocessing of Xrays into longer wavelengths for an individual system, and more sophisticated modelling is needed to explain the origin of the optical/UV emission during outbursts (Vrtilek et al. 1990;Gierliński et al. 2009).
On the other hand, these studies have been traditionally performed for LMXBs, since for high-mass X-ray binaries (HMXBs) the emission from the optical companion is expected to dominate over reprocessing effects at optical and ultra-violet (UV) wavelengths.Indeed, such bright optical outbursts taking place simultaneously with the X-ray outbursts are not common in BeXBs.As far as we know, the only system displaying very bright optical outbursts during active X-ray epochs that has been reported, is the Be-ray binary 1A 0538-66 (this source also reached L X ∼ 10 39 erg/s), which can reach amplitudes of the variability of ∼2 mag in the V-band (Charles et al. 1983).Several mechanisms have been proposed to explain the extremely bright optical outbursts of this system, including irradiation inside the binary system, although this system showed many other peculiarities with respect to other BeXBs, which complicated the interpretation of the observations (Charles et al. 1983).Ducci et al. (2019) performed a very detailed analysis of the optical emission of this system, considering both irradiation of the accretion disk and heating of the surface of the Be star, but concluded that in the case of 1A 0538-66, the optical emission is mainly powered by reprocessing of X-rays in the envelope surrounding the binary system.
In order to constrain the origin of the optical/UV emission, there are several methods that can be used: -Compare the light curves in the different optical/UV bands with the X-ray light curve.
-Study the correlation between the measured optical/UV and the X-ray fluxes.-Measure the timescales of the exponential decay of the outburst for the light curves in each band.-Develope physical models in order to reproduce the observations.
In this paper, we have carried out a detailed analysis of the optical and UV emission from Swift J0243.6+6124, in correlation with the X-ray fluxes during the 2017-2018 giant outburst, performing the calculations mentioned above.In Sect.2, we describe the observations used on this analysis.In Sect.3, we present the calculations and data analysis.In Sect.4, we discuss the possible origins of this emission and describe the physical model we propose to reproduce the optical/UV emission.Finally, in Sect.5, we summarize the main results and present the conclusions of the paper.

Optical/UV photometry
The optical V-Johnson band light curve contains data from several observatories (see Reig et al. 2020 for a detailed description of the optical data).We included the light curve from the ASAS-SN Variable Stars Database (Shappee et al. 2014;Jayasinghe et al. 2019), correcting the data from the contamination due to the presence of a close star in the photometric extraction aperture.Photometric observations from Skinakas and Aras de los Olmos observatories are also included.There are also some observations from the American Association of Variable Star Astronomers (AAVSO; Watson et al. 2012).Additionally, we have also included the photometric observations in the v filter from the Ultra-Violet/Optical Telescope (UVOT, Roming et al. 2005) on board the Neil Gehrels Swift Observatory (hereafter Swift; Gehrels et al. 2004) instrument.To convert the UVOT-v filter measurements to the Johnson V-band we have applied a filter conversion correction.This correction has been calculated by comparing the synthetic photometry of the theoretical spectrum of an O9.5V star reddened with E(B-V)=1.20 mag (Reig et al. 2020) in both filters.
The long-term optical light curve of Swift J0243.6+6124 in the V-Johnson band, together with the Swift/BAT (15-50 keV) light curve, is shown in Fig. 1.
Multi-colour photometric observations from the other UVOT filters have also been analysed in this work.The UV and optical light curves of Swift J0243.6+6124 were obtained after checking the images using the uvotsource ftool to do aperture photometry.Data where the source fell on a low-sensitivity patch have not been included.

Swift/XRT and BAT observations
The Swift observatory monitored the outburst of Swift J0243.6+6124 in detail.In order to characterise the X-ray evolution throughout the outburst, we used the data from both the X-ray Telescope (XRT; Burrows et al. 2005) and the Burst Alert Telescope (BAT; Barthelmy et al. 2005).The XRT performed pointed observations of the outburst at a variable but high (up to daily) cadence, changing between photon counting (PC) and window timing (WT) modes throughout the outburst and its rebrightening events dependent on the target's brightness.A large fraction of the observations was taken in automatic mode, accumulating the majority of data in the suitable mode for the source flux levels.We extracted XRT X-ray target and background spectra, as well as response and ancillary files, using the Swift online data reduction pipeline (Evans et al. 2009, www.swift.ac.uk/ user_objects/) for all observations of the 2017/2018 outburst.This pipeline automatically accounts for the pile-up during the brightest phases of the outburst.For Swift/BAT, we downloaded the daily source count rates from the Hard X-ray Transient Monitor webpage (Krimm et al. 2013, https://swift.gsfc.nasa.gov/results/transients/).We note that, compared to the BAT light curves as reported during the outburst, any anomalously low count rates due to saturation during the outburst peak have been corrected in the current version of the monitor webpage.

Data analysis
In this section, we present the methods that we have used to convert from observed photometry to monochromatic fluxes, in order to build the Spectral Energy Distributions (SEDs) corresponding to each time.We also present a new estimation of the interstellar reddening.After that, we describe the method to isolate the extra emission intrinsic to the outburst.The method used to estimate the bolometric X-ray fluxes is also described.Moreover, in order to compare with previous multi-wavelength analysis (see Sect 1), we calculate the timescale of the exponential decay in the different bands, and the F X /F opt/UV correlation indexes.

Transformation from observed magnitudes to monochromatic fluxes
To convert the observations in the V-Johnson band to monochromatic fluxes, we used the standard zero point method.The UVOT magnitudes are provided in the AB-photometric system.
In order to convert them to monochromatic fluxes and be able to use these measurements to build the SEDs, we calculated the corrections to the filter response depending on the spectral shape of the observed object, as described in Brown et al. (2016).As these authors explain, the monochromatic fluxes provided by the pipeline are not always correct, since these values are obtained using the zero points corresponding to Vega and this difference can be crucial at the UV wavelengths.We calculated the flux correction (FC) for each filter, considering the theoretical spectrum of an O9.5V star, and reddening of E(B-V)=1.20 mag (see Sect. 3.2).The visualization of the differences between considering Vega or the theoretical reddened spectrum of Swift J0243.6+6124 can be observed in Fig. 2. The largest differences are found for the UVOT-uw1 and UVOT-uw2 filters.

Estimation of the interstellar reddening
In order to identify the contribution of each emission mechanism to the total optical and UV observed light, correctly determining the reddening is needed.
In Paper I, we obtained the spectral type of the optical companion to be an O9.5 V star, and we provided two independent measurements of the extinction: from optical photometric observations we obtained E(B-V)=1.24±0.02mag, while when we considered the strength of the diffuse interstellar bands (DIBs) we got E(B-V)=1.1±0.2 mag.As we mentioned in Paper I, the determination of the interstellar extinction in Be stars (and in BeXBs) can be affected by the presence of the circumstellar disk around the Be star, which leads to an excess in the IR and redder optical wavelengths.We present here a new determination of the interstellar extinction, obtained by fitting the observed Spectral Energy Distribution (SED) points, including also the UV data from Swift/UVOT, with synthetic photometry calculated from a grid of reddened theoretical spectra.
In order to minimize the contribution from the circumstellar disk and the emission related to the X-ray outburst, we selected the faintest observations from the long-term light curves, obtained in the least active X-ray epoch (close to the epoch around MJD 58530, see Fig. 1).In this case, we would only expect to observe the contribution of the Be star and the circumstellar disk.From the long-term variability of the optical light curves, and from the values of the equivalent width (EW) of the Hα line, we expect that the contribution from the Be circumstellar disk is minimum at that epoch (Reig et al. 2020;Liu et al. 2022b).For the SED fit, we built a grid of values of E(B-V) from 0.9 to 1.4 with step 0.01, and used the theoretical spectra of a O9.5 V star of T e f f =32000 k, log g = 4.0 and Z=Z⊙ from the Castelli & Kurucz 2004 Stellar Atmosphere Models (Castelli & Kurucz 2003), assuming the extinction law from Fitzpatrick (1999).We considered a distance of 5.2 kpc and obtained synthetic photometry to compare with the observations.We got the best fit for E(B-V)=1.20±0.01mag, consistent with the determination in Paper I. Similar results are obtained if the SED fit is performed only to the observed fluxes in the U+UV filters (excluding the B, and V filters, which would be more affected by the contribution of the circumstellar disk).The results of the SED fit are shown in  in agreement with the evolution of the EW(Hα) line (Reig et al. 2020;Liu et al. 2022b).In order to have an estimation of the emission from the Be star and the circumstellar disk in the V-Johnson band, we performed a cubic spline fit to the points outside the outburst (we did not consider the points between MJD 57850 and 58300 MJD) on the long-term V-Johnson light curve (see Fig. 5).For the other optical and UV bands for which the temporal coverage was poorer, we just performed a linear fit to those points out of the outburst between MJD 58200-58300 (epoch when it is considered that the giant outburst ended) and MJD 58500-58600 (epoch when the minimum fluxes were measured at all the wavelengths).Then we substracted this linear fit to the observed fluxes along the outburst.The long-term light curves of the fluxes in the six UV/optical filters are shown in Fig. 4. The estimation of the Be+disk emission at each wavelength is plotted with a gray dashed line.The variations due to changes in the circumstellar disk are larger for redder wavelengths (as can be seen from the slopes of the fits in Fig. 4).

Isolation of the extra optical/UV emission due to the outburst
We used these residuals (after substracting the Be+disk component) to build the daily optical/UV SEDs that trace the emission exclusively due to the outburst.These SEDs will be shown and studied in Sect.4.2.2.

Estimation of the X-ray fluxes
We fitted each XRT spectrum with a combined absorbed power law plus black body model, scripting the automatic fit for each spectrum in xspec (Arnaud 1996).We assumed interstellar abundances from Wilms et al. (2000) and cross-sections from Verner et al. (1996).We used the tbnew, pegpwrlw, and bbodyrad models of xspec.In all fits, the neutral hydrogen absorption column density was fixed to N H = 0.825 × 10 22 cm −2 .We have calculated this quantity considering a variety of relations from the literature.We applied the linear relations between N H and A V given by different authors.We calculated the value of the optical extinction A V = R V *E(B-V)=3.7 mag, from the value of E(B-V)= 1.20 mag obtained in Sect.3.2 and considering R V =3.1.From the relation obtained by Güver & Özel (2009), we got N H = 0.822 × 10 22 cm −2 .A value of N H = 0.818 × 10 22 cm −2 is obtained when using the relation given by Watson (2011).A slightly higher value of N H = 1.07 × 10 22 cm −2 is obtained if we use the relation given by Foight et al. (2016).On the other hand, the relation given by Willingale et al. (2013) from the value of N HI = 0.726 × 10 22 cm −2 (got with nh of FTOOLS) would provide an estimation of N H = 0.87 × 10 22 cm −2 .While more complex spectral models have been used in modelling the X-ray spectra of other observatories (see e.g.Zhang et al. 2019;Jaisawal et al. 2019), this approach suffices for an accurate flux determination in the Swift band.From the Swift/XRT spectral fit, we obtained the fluxes in 0.5-2 keV and 2-10 keV bands that will be used to calculate the correlation indexes in Sect.3.6.
To obtain bolometric fluxes from the Swift monitoring, however, the spectral cutoff at higher energies needs to be corrected.Therefore, we combined the above Swift/XRT fits with the measured 15-50 keV count rates from the Swift/BAT hard X-ray transient monitor.As the spectral shape changes during the outburst, we cannot simply apply a single, constant conversion factor between the BAT count rate and flux.Instead, we took the best-fit XRT spectral model for each observing date and then added a high-energy exponential cutoff (highecut model in xspec).So that the total model has the functional form of: Here K is the normalization of the blackbody component, A is the normalization of the power law, and α is the photon index of the power law.The relevant parameters for the high energy cutoff are E c , the cutoff energy, and E f , the folding energy, both in keV.At energies below E c , the model is equivalent to the absorbed power law plus blackbody spectrum used for the XRT fits.We assumed that the exponential cutoff only starts at 8 keV (i.e., we fixed E c = 8 keV), with a variable characteristic rollover or folding energy E F .This way, any change in the folding energy that is necessary to describe the BAT spectrum, does not influence the XRT data, as they were not used above 8 keV.At each attempted folding energy, we calculated the corresponding BAT count rate and compared this value to the measured value on that day.Storing the folding energy yielding the best agreement between cutoff-model and observations, we finally calculated the bolometric (0.1 − 500 keV) X-ray flux using the modified spectral model.From this fit, we also obtained the fluxes in the 15-50 keV band that will be used in Sect.3.6.
The light curves of the residual optical emission in the V-Johnson band obtained from the calculations described in Sect.3.3, along with the X-ray light curve of the Swift J0243.6+6124outburst, converted to luminosities considering a distance of d=5.2 kpc, are shown in Fig. 6.A clear correlation between the optical and X-ray fluxes is observed, although this correlation differs between the rise and the decay of the outburst.

Timescales of the light curves decays
As mentioned in Sect. 1, when reprocessing dominates the optical emission during X-ray outbursts of X-ray binaries, the exponential decay timescales of the optical and IR light curves should be 2-4 times longer than the X-ray value (King & Ritter 1998;Rykoff et al. 2007).An exponential decay can be observed for Swift J0243.6+6124 in Fig. 6.
We have calculated the decay of the outburst in our Xray, UVOT-uw2, UVOT-um2,UVOT-uw1, and V-Johnson light curves, by fitting the formula F = F 0 e −t/τ to our observations, and the results of these fits are shown in Table 1.We could not fit the UVOT-u and UVOT-b light curves due to the poor temporal coverage of the outburst decay that did not allow to correctly perform the fit.
Table 1.Timescales of the decay of the giant outburst for the X-ray, UVOT-uw2,UVOT-um2,UVOT-uw1, and V-Johnson light curves, as a result of fitting the formula F = F 0 e −t/τ to the observed light curves.
The decay timescales of the optical (V-Johnson band) and UVOT-uw1 light curves are ∼4 times the decay timescale of the X-ray light curve.However, the decay timescales of the UVOT-uw2 and UVOT-uw1 light curves are ∼2 times the decay timescale of the X-ray light curve.

Flux-flux correlations
As mentioned in Sect. 1, a power law correlation F V ∼ F 0.5 X have usually been interpreted as reprocessing of X-rays at longer wavelengths in X-ray binaries (van Paradijs & McClintock 1994;Russell et al. 2006).
We have calculated the indexes of the correlation between the fluxes in each optical/UV filter and simultaneous X-ray fluxes, considering different X-ray energy bands.The correlations between the optical/UV logarithmic fluxes in the different filters and the bolometric logarithmic X-ray fluxes (0.1-500 keV) are represented in Fig. 7.For the fluxes in the redder wavelengths (V-Johnson until UVOT-uw1 filters) a different correlation is observed during the rise than during the decay of the outburst.For the V-Johnson and UVOT-uw1 filters, two fits (for the rise and for the decay) have been performed.Although a different trend is also observed during the rise than during the decay for the observations in the UVOT-b and UVOT-u filters, the temporal coverage for the decay of the outburst was not good enough to perform a reliable fit, and only the fits during the rise of the outburst are provided.For the filters at shorter wavelengths (UVOT-um2 and UVOT-uw2) this effect is not observed and a global fit provides better results.The results of this correlation analysis are provided in Table 2.
Appart from the different trends between the rise and decay for the redder fluxes, we have also observed that the correlation indexes depend on the energy band used to estimate the X-ray flux, increasing when harder X-rays are used.This effect is observed for all the optical/UV filters (see Table 2).

Results and discussion
In this section, we will discuss first the possible mechanisms that can explain the observed UV/optical fluxes, considering the results obtained in Sect.3. Then we will model the most relevant effects in order to reproduce the evolution of the SED along the outburst in correlation with the X-ray emission.We will discuss the possible interpretations of the results of the SED fits.Finally, we will provide some possible explanations to the optical/UV humps observed during the decay of the outburst.

Phenomenological constraints on the optical/UV origin
As we mentioned in Sect 1, there are several mechanisms that can contribute to the optical/UV emission and display similar amplitude variability in X-ray binaries.
Firstly, we will consider the mechanisms inherent to the Be star and/or the circumstellar disk.Long-term variability, either due to formation and dissipation of the Be circumstellar disk and/or to changes in the inclination or due to warping and/or precession of this circumstellar disk, typically produce optical variability with amplitudes of tenths of magnitudes (Yan et al. 2012;Moritani et al. 2013;Alfonso-Garzón et al. 2017), but the timescales of this circumstellar disk variability are longer than the duration of the X-ray outbursts of BeXBs.Indeed, X-ray outbursts usually take place when the circumstellar disk reaches a critical size, and/or the inclination of the circumstellar disk favours the accretion onto the NS.On the other hand, optical outbursts due to mass ejections from the Be star with durations of several days, and displaying amplitudes of tenths of magnitudes are commonly observed in classical Be stars (Rivinius et al. 2013) and in BeXBs (Yan et al. 2012;Alfonso-Garzón et al. 2017).However, when ejections of matter from the donor star take place, a delay between the X-ray emission and the optical outburst would be expected.Indeed, they can cause the replenishment of the Be circumstellar disks, and these changes in the Be disk could trigger an X-ray outburst (Yan et al. 2012).Then, these mechanisms can not explain the tight correlation and strict simultaneity between the optical/UV fluxes and X-ray fluxes that we observe during the giant outburst of Swift J0243.6+6124.
We also considered that the optical emission could be originated in the jet, given the fact that this was the first Be X-ray binary for which signatures from a jet has been found (van den Eijnden et al. 2018).We have estimated the contribution of the optical emission considering the radio spectral shape, which is steep at the peak of the outburst (α<0).If we take the brightest observed radio flux at 58073 MJD, which was 92 µJy at 6 GHz, an we use the spectral index measured by van den Eijnden et al. ( 2018) α=-0.64,we obtain an extrapolated optical flux of F V ∼ 6.1 × 10 −20 erg s −1 cm −2 Å −1 .If we do the same calculation for the observations at 58090 MJD, when the slope Article number, page 7 of 16 A&A proofs: manuscript no.SwiftJ0243_outburst Table 2. Slopes of the correlations between the optical/UV fluxes and the X-ray fluxes in different bands, F opt/UV ∝ F β X , during the outburst.was flatter, α = 0.08, we obtain an optical contribution from the jet of F V ∼ 6.2 × 10 −17 erg s −1 cm −2 Å −1 , still very small compared to the observed variability intrinsic to the outburst of ∆F V ∼ 2.3 × 10 −13 erg s −1 cm −2 Å −1 (see Fig. 7).We can then rule out the emission coming from the jet as the main origin of the optical/UV emission.Moreover, there were also observations after the optical outburst with similar radio fluxes when the optical fluxes were fainter, so that we can conclude that the optical and radio emission are not correlated in this system.We want to point out that the BeXB radio jets are really faint: ∼6200 times fainter than in black hole X-ray binaries (see van den Eijnden et al. 2022).That could explain why you may possibly have jet contributions to the OIR emission in BH LMXBs, but not for BeXBs.
There are of course, other mechanisms that could lead to optical/UV emission directly correlated with the X-ray emission and have to be considered in order to explain the optical and UV fluxes observed during this outburst.These mechanisms are: a) thermal emission from the viscously heated accretion disk around the NS, b) reprocessing of the hard X-rays in the accretion disk, c) X-ray heating of the surface of the Be star.Given the tight correlation observed between the X-ray and optical/UV fluxes, some of these X-ray related mechanisms could be significantly contributing to the optical and UV emission observed during the giant X-ray outburst of Swift J0243.6+6124.The UV and optical timescales calculated in Sect.3.5 are in the range 2-4 times the X-ray decay timescale (see Table 1), so according to King & Ritter (1998), these results favour the X-ray reprocessing as the main mechanism contributing to the UV and optical variability observed during the outburst.The differences found between the timescales of UVOT-uw2 and UVOT-um2 and those obtained for UVOT-uw1, and V-Johnson (the last ones are around twice longer than the other ones), could be indicating that there are differences between the origin of the emission on the different bands.The correlation indexes obtained from the flux-flux correlation analysis (F V ∝ F β X ) are shown in Table 2. Some of the obtained indexes, would be in agreement with the correlation index typical of X-ray reprocessing proposed by van Paradijs & McClintock (1994), β=0.5.However, the obtained correlation indexes vary between the different optical/UV filters and when different X-ray bands are used.Moreover, for the redder observations, these indexes are larger for the rise than for the decay of the outburst.The dependence on the energy band could be related to the changes in the X-ray spectra, and the variable contribution from the reflected emission during the Super-Eddington state (Bykov et al. 2022).Also the fact that the albedo is expected to vary with the energy of the X-ray impinging radiation could affect to the correlations (harder X-rays penetrate more in the disk).The differences between the rise and decay could be related to the presence of a more extended accretion disk during the decay than during the rise, or if the disk gets flared along the outburst.Moreover, the source underwent different accretion regimes, from sub-to super-critical accretion, and to a third regime with a radiation-pressure dominated accretion disk (Doroshenko et al. 2020), and this could also affect to the correlations.As mentioned by Gierliński et al. (2009), these simple flux-flux correlations alone can be inaccurate to explain in detail the origin of the optical/UV emission for individual LMXBs.For HMXBs, the analysis can be even more complicated.For this reason, we tried a different approach involving a more complex modelling of the emission mechanisms that will be presented in Sect.4.2.
We also wanted to give a possible explanation to the singularity of this optical/UV outburst.The unprecedent high X-ray luminosities for a BeXB observed during this giant outburst, being the first Ultraluminous X-ray Pulsar (ULXP) identified in the Galaxy, could have led to the exceptionally bright outburst at optical and UV wavelengths.X-ray outbursts of BeXBs typically have luminosities of L X =10 36 -10 37 erg/s.For the orbital parameters of Swift J0243.6+6124, and considering the intrinsic brightness of the optical companion, an X-ray outburst of L X =10 37 erg/s would produce an observed optical outburst with amplitude ∆V=0.06 mag, and an X-ray outburst of L X =10 36 erg/s would lead to an optical outburst with ∆V=0.02mag.These values are close to typical uncertainties in optical photometric measurements in many cases, and are much smaller than the amplitude we observe in our optical outburst.On the other hand, the orbital period of Swift J0243.6+6124(P orb =27.587 days), and hence the distance between the NS and the Be star, is shorter than for other BeXBs that have underwent very bright giant outbursts, such as 1A 0535+262, which reached L X ∼ 1 ×10 38 erg s −1 (Reig et al. 2022), and whose orbital period is P orb =111.1 d (Finger et al. 2006).Then, we can explain the exceptionality of this optical/UV outbursts as a combination of the extreme X-ray luminosities, together with a semimajor axis value shorter than for other BeXB systems.

SED modelling through the outburst
Since the effects directly related to the X-ray emission seem to be the most plausible mechanisms driving the optical/UV emission observed during the giant ourburst of Swift J0243.6+6124,we have created a physical model considering the following mechanisms: emission from gravitational energy released through viscosity, emission coming from the X-ray irradiation of the accretion disk, and emission due to the heating of the surface of the companion star.We have modelled these effects and used the calculated fluxes in each optical/UV filter to build a theoretical SED to be compared with the observed optical/UV SEDs.

Description of the model
In all the calculations, we have used a distance of d = 5.2 kpc to estimate the X-ray incident luminosities (which is an input for the three models), and to convert the calculated optical/UV luminosities to fluxes for comparison with the observed SEDs.
We have used the most recent orbital parameters provided in the GBM Accreting Pulsar Histories web accesible at https://gammaray.nsstc.nasa.gov/gbm/science/pulsars/lightcurves/swiftj0243.html of the GBM Accreting Pulsars Program (GAPP, Malacaria et al. 2020).We then used an orbital period P orb = 27.698899days, a projected semimajor axis a sini=115.531light-sec, an eccentricity e=0.1029, and a binary epoch T π/2 =2458116.09700JED.On the other hand, we considered an inclination of the orbit of i=30 • (Paper I).
For each day of the outburst (whenever we had some optical/UV observation), we used the X-ray luminosity extracted from linear interpolation of the X-ray bolometric fluxes in the light curve generated from the method described in Sect.3.4, to be used as an input of the different parts of the model.The XRT spectra (purple line) and the models used to fit them for two representative days of the outburst, one close to the peak of the outburst (58062 MJD, up) and another one during the decay of the outburst (58094 MJD, bottom), are plotted in Fig. 8.The black- For the two mechanisms involving the accretion disk, we performed two different fits.In the first one, we left the outer radius R out as a free parameter, using a grid with step 0.01a, and taking values from R out =0.01a to R out =0.20a (typical predicted sizes of an accretion disk according to theoretical models are R out = 0.04 − 0.10a, see for example Hayasaki & Okazaki (2004, 2006))), being a the semi-major axis of the orbit.We have set an upper limit of R out based on the size of the Roche lobe.For a mass of the donor star of M * =18 M⊙ and considering the mass of the NS M NS =1.4 M⊙, and using the approximation by Eggleton (1983), we got a value of the Roche lobe of R L =0.20a.For the second fit, we fixed R out = 0.20a and fit a factor dependent on other effects that will be explained later.For each accretion rate, we approximated the inner radius of the accretion disk as the magnetospheric radius, calculated considering a magnetic field of B∼6.6*10 12 G (value obtained from the estimation by Wilson-Hodge et al. (2018), B∼ 1 × 10 13 (d/7kpc) G).To estimate r m we used the following equation (Fürst et al. 2017): (1) The accretion rate ( Ṁ) has been obtained from the X-ray bolometric luminosity using the formula: and considering typical values of a NS of M = 1.4 M ⊙ and R = 10 6 cm.
We describe below the details of the specific models we considered for each of these three mechanisms.These components of the model used to fit the optical/UV observations (red points), are plotted in Fig. 8.
-Emission from the viscously heated accretion disk: Intrinsic thermal emission from the viscously heated outer accretion disk can contribute to the X-ray, UV and optical emission of X-ray binaries (Shakura & Sunyaev 1973;Frank et al. 2002).We modelled the contribution of a viscously heated accretion disk, using the approximation of an optically thick disk (Shakura & Sunyaev 1973), using the formula: (3) We calculated the temperature for each radius between R in and R out , and integrated the Planck law over each area element of the disk to estimate the emitted spectrum.
The modelled emission due to this effect is plotted with green dash-dotted lines in Fig. 8.As can be seen, due to the effect of the magnetic field on the inner radius, this emission does not contribute to the soft X-rays in the energy range covered by the XRT observations (purple line).
-X-ray irradiation of the accretion disk: We have modelled the effect of the X-ray irradiation of the accretion disk for each incident L X along the outburst, for the grid of values of R out described above.We have modelled the emission from an irradiated disk using the calculations in Vrtilek et al. (1990), which give this relation for T(R): where f 1 is a factor that depends on the details of the vertical structure of the disk and has been fixed to 1, f 2 is the absorbed fraction of the impinging radiation and has been fixed to 0.5, and f 3 takes into account the possible anisotropy of the radiation and has been fixed to 1.The modelled emission due to this effect is represented with brown dashed lines in Fig. 8.
-X-ray heating of the surface of the companion star: The irradiation of the surface of the the Be star by the Xray photons emitted close to the NS can contribute significantly to the optical and UV emission (van Paradijs & Mc-Clintock 1995;Charles & Coe 2006).We have calculated the extra emission in each band by estimating the difference between the emission from an unirradiated O9.5 V star and that from an irradiated star, considering a radius of the donor star of R * = 8R ⊙ , and an initial effective temperature of T eff =32000 K.
We considered that the X-rays emitted by the NS are absorbed by the atmosphere of the optical companion and then thermally re-radiated at lower energies.We have followed a similar modelling to the one describe in Ducci et al. (2019) to estimate the contribution of the heated surface of the Be star on A0538-66.In this model, the X-rays emitted by the Article number, page 10 of 16 central source are absorbed by the atmosphere of the companion star and heat its surface.For each surface element of the star, the local effective temperature increases as: with T 4 0,e the effective temperature of an unirradiated surface element, η is the albedo (fraction of reflected X-rays) that we have fixed to 0.5, L X is the X-ray luminosity of the source, ϕ is the angle between the direction to the X-ray source and the normal of the surface element, σ is the Stefan-Boltzmann constant, and R is the distance from the surface element to the X-ray source.We calculated this distance for each phase using the latest orbital parameters.We have calculated the contribution of this effect for each Xray luminosity along the outburst.To illustrate this effect, the modelled heating of the surface of the optical companion of Swift J0243.6+6124 at the peak of the outburst (MJD 58064) is shown in Fig. 9.As can be seen, the center of the facing star reaches an extra heating of ∆T > 12000 K.We have tried to simplify this contribution of the model as much as possible and we are aware that there are some considerations that we have not taken into account.We want to point out that, since the temperatures would be always larger than the effective temperature of the star, the emitted spectrum would always peak at wavelengths shorter than 905 Å (λ peak for T e f f =32000 K from Wien's law).This corresponds to energies beyond those in the range of our optical/UV observations, which all fall on the Rayleigh-Jeans tail of this emission (see Fig. 8, the emission reprocessed by the Be star is plotted as a pink dotted line).This implies that systematic uncertainties on the details of the reprocessing process would mainly have a small impact on the normalization only (increasing or decreasing the modelled continuum emission), but not on the shape of the spectrum itself.
One example of such an uncertainty is the exact value of the albedo, which affects the range of temperatures reached by the heated stellar atmosphere (see Eq. 5).We have made several tests modifying this value and have concluded that, compared to the model represented in Fig. 8, these modifications would cause only small changes in the normalization, but not in the spectral shape at the observed wavelengths.Another uncertainty is the role of the stellar rotation: Eq. 5 implicitly assumes that the stellar atmosphere is able to reach thermal equilibrium with the incoming X-ray flux on timescales significantly shorter than the rotational timescale of the star.In Swift J0243.6+6124, that timescale is of the order τ rot ∼ 1 day, if we assume the values of the projected rotational velocity vsini=210 km/s and inclination i=30 0 from Reig et al. (2020).From the correlation observed between the optical and X-ray observed fluxes, we can conclude that the cooling/heating timescale should be shorter than the rotational period.But if the timescale was of the order or longer than the rotational period, the X-rays would produce a temperature profile with constant-temperature rings at each latitude, and with peak temperatures slightly lower than when no rotation is considered (due to the redistribution of the incoming energy).In any case, both faces of the star would be heated, although they would reach maximum peak temperatures slightly lower (e.g.∼40700 K vs. 44700 K in the equator of the Be star at the peak of the outburst), and the modelled emission would decrease by a factor of ∼0.8 in the spectral range of our observations.On the other hand, we are not considering that some fraction of the donor star could be obscured by the circumstellar disk, although the height scale is expected to be small (Okazaki et al. 2002).Modifying the fraction of the surface of the Be star that is heated would also lead to a variation in the normalization of the blackbody emission.Detailed calculations of these timescales, as well as extra effects such as differential rotation and convection that further alter the temperature profile, or modelling the occultation by the circumstellar disk would introduce a large set of new assumptions, compared to the limited number of fitted datapoints.We therefore leave the consideration of these effects to future modelling work with more complete datasets.
We also want to mention that we have not included in our modelling the effect of the irradiation of the Be circumstellar disk.It should be taken into account that there are some important differences between circumstellar disks in BeXBs and Be disks in classical isolated Be stars.The correlation between the orbital period in BeXBs and the H α equivalent width has been interpreted as evidence of the truncation of the Be disk (Reig 2011).These disks can be truncated at different resonant radii, where the angular frequency of the disk and that of the orbit is an integer (Monageng et al. 2017), and can reach the size of the Roche lobe during X-ray outbursts (which is R L ∼10 R * for the donor star in this system).However, due to the resonant torque the outer disc density drops much more rapidly in Be X-ray binaries than in isolated Be stars (Okazaki et al. 2002).
In our modelling, we assumed that the Be disk was coplanar to the orbit during the outburst, which is the simplest configuration to explain the occurrence of an X-ray outburst in a BeXB system (Okazaki et al. 2002).We are aware that some proposed scenarios to explain the occurrence of giant outbursts involve different configurations, including precession of warped Article number, page 11 of 16 Be disks (Martin et al. 2011) and eccentric and tilted disks due to Kozai-Lidov oscillations (Martin et al. 2014;Martin & Franchini 2019).Indeed, for some BeXBs, signatures of tidal warped circumstellar disks have been observed during giant X-ray outbursts (Moritani et al. 2013;Reig & Blinov 2018).However, no signs of warping and/or eccentricity have been found in the profiles of the optical emission lines during the giant outburst of Swift J0243.6+6124(Reig et al. 2020;Liu et al. 2022b), although we can not rule out the presence of a tilted disk with the tilting axis almost perpendicular to the line of sight (in this case the emission lines would also be single-peaked).
Under the considered assumption that the Be disk is coplanar to the orbit, the X-rays would not impinge on the Be disk, since the angles that the X-rays form with the normal vector of the Be disk surface would be close to 90 degrees and then they would not heat the Be disk surface.Moreover, a significant fraction of the X-rays inciding with such small angles w.r.t. the orbital plane could be firstly intercepted by the outer parts of the flared accretion disk.Then, under this assumption, only the outer, low-density, and colder material of the disk could be heated.The emission of this colder material, if irradiated, would probably peak at redder wavelengths.We consider that correctly modelling the temperature and density gradient and the geometry of the Be disk, and the effect of the X-ray irradiation of this disk would require many more assumptions, and that performing such a detailed analysis with the dataset we present in this work is out of the scope of this paper.However, we want to point out that not including this contribution could lead to an overestimation of the outer radius of the accretion disk (see Sect. 4.2.2),since similar temperatures could be reached in the outer regions of both disks.

Optical/UV SEDs fits
The optical/UV SEDs built with the available observations along the outburst, together with the best SED fit for each day are plotted in Fig. 10.For those days for which we had optical measurements in the V band, we calculated the value of the outer radius of the accretion disk, R out , yielding the best fit of the SED using the minimization of chi-square method (see Fig. 11).The black solid lines represent the sum of the three modelled components.The different mechanisms are plotted with different colours and line styles: the pink dotted lines represent the modelled contribution by the extra emission from heating the surface of the companion star, and the brown dashed lines correspond to the emission modelled from an irradiated accretion disk with the value of R out providing the best fit.The emission from the viscously heated accretion disk is too faint to appear represented in the plot.Since it does not contribute significantly to the optical and UV fluxes and can be ruled out from the discussion.
As can be seen in Fig. 10, the heating of the surface of the companion star (blue dotted line) contributes significantly to the optical and UV observed fluxes along the whole outburst.The SED fits get improved with the addition of an irradiated disk from MJD 58047 to 58120.Before MJD 58047, and from MJD 58141 onwards, the heating of the companion star is enough to explain the observed optical/UV emission.Although from MJD 58141 onwards, there were not enough optical observations to fit the redder part of the SED (more affected by irradiation of the accretion disk at these X-ray luminosities).The irradiation of the accretion disk contributes significantly to the UV emission close to the peak of the outburst (from MJD 58062 to 58096 aprox.).
As mentioned above, to reproduce the observed fluxes, we performed two different fits.First, we left the outer radius of the accretion disk (R out ) as a free parameter of the fit (see Fig. 11).From the results of the fits, an increasing R out along the outburst is obtained, reaching a value of R out ∼0.18a at 58097 MJD.After that, a larger value of R out , reaching the upper limit given by the Roche Lobe R out ∼ R L = 0.20a, is required to fit the observed optical fluxes from 58100 to 58120 MJD (see Fig. 11).This could be understood as an increase of the size of R out and could be interpreted as a progressive increase of the size of the accretion disk or a dissipation towards the outer parts of the disk.We want to stress that for this first attempt, we fixed all the geometrical factors and the fraction of reprocessed emission along the whole outburst, and we only considered R out as a free parameter in our modelling.This increase of the emission contributing to the redder wavelengths along the outburst is in agreement with the differences observed between the rise and decay in the fluxflux correlations for the redder wavelengths (Sect. 3.6).
Such an increase on the redder fluxes could also be reproduced by modifying the factors in the irradiation disk model which depend on the vertical scale of the disk and the effects of possible anisotropy of the radiation, or the fraction of reprocessed emission due to changes in the ionization or the spectra shape (see Eq. 4).Such effects seem to have occurred during this outburst.Doroshenko et al. (2020) proposed that during this giant outburst, the system underwent transitions between three states: sub-critical gas pressure dominated disk (GPD), supercritical GPD (dashed blue line in Fig. 6), and super-critical radiation pressure dominated disk (RPD).The times proposed for the transitions between these states are marked as dashed and solid pink lines respectively in Fig. 6.These changes in the geometry of the accretion disk should have an effect on the fraction of reprocessed X-ray emission at longer wavelengths.Moreover, the change from pencil to fan beam regime proposed by Liu et al. (2022a) also should have a direct impact on the fraction of Xray flux impinging on the accretion disk due to geometrical effects.It could also happen that the accretion disk gets flared, and this would result in a stronger heating of the outer parts of the disk and a larger contribution to the redder fluxes.Moreover, the X-ray radiation was harder during the decay than during the rise of the outburst (for similar X-ray fluxes), and the albedo is smaller for harder X-rays.A smaller albedo would produce higher optical/UV emission.To quantify these effects, we have also performed the SED fits fixing R out =0.20a, and leaving as free parameter a combined factor k = f 2 f 3 f 1 , which includes the three factors present in Eq. 5.The results of these fits are shown in Fig. 12.As expected, and increase of this value is obtained along the outburst.
On the other hand, it is believed that the accretion disk displayed outflows during the brightest part of the outburst (van den Eijnden et al. 2019b), and an extra contribution to the optical continuum emission could be originated in the base of this disk wind (Koljonen et al. 2023).Finally, the timescales of the cooling of the accretion disk and the donor star could have also an effect on the observed evolution along the outburst.

Optical hump at around MJD 58095-58112 and end of the outburst
If we look in detail to the optical/UV light curves of Swift J0243.6+6124, a hump from around MJD 58095 to 58112 is observed during the decay of the outburst at all the optical wavelengths and also in the UVOT-uw1 filter (see Fig. 4).
Article number, page 12 of 16    In the previous section, we mentioned that the observed increase on the contribution of the irradiated accretion disk along the outburst could be related with changes in the geometry.Liu et al. (2022a) studied the variability of the X-ray pulse profiles along the outburst and identified a high-double peaked pulse profile present from MJD 58048 to 58090, and called it the accretion column regime (see blue solid lines in Fig. 6).Doroshenko et al. (2020) proposed that a transition between the GPD to RPD disk and viceversa took place at around MJD 58045 and 58098 (pink solid lines in Fig. 6).During this epoch, we observed a progressive increase of the R out and k values obtained from the SED fits (see Fig. 11 and Fig. 12).After MJD 58097, a sudden change in the fitted values is observed.According to Liu et al. (2022a), the pulse-profile gradually changed from double-to single-peaked from MJD 58091 to 58105 when the pulse profile was single peaked again (dashed blue line in Fig. 6).This transition epoch would be roughly coincident with the epoch when we observe the optical/UV humps and the change in the fitted R out and k is observed.On the other hand, if a transition from RPD to GPD took place at around MJD 58098, it should be related with the plateau of the optical and UV fluxes.A sudden change in the geometry of the inner flow would be also in agreement with the change on the fitted values of the combined factor k at around MJD 58098 (see Fig. 12).If the optically thick curtain vanishes, there should be an increase of the X-ray irradiation impigning the accretion disk that could produce the observed humps from MJD 58095 to 58115 and the sudden increase on the fitted values of k.
From MJD 58105 to 58139, the source transited from the single-peak pulse profile to the low double-peaked profile, that has been interpreted as the transition from accretion column to the low pencil beam regime (Liu et al. 2022a).The impinging X-ray radiation onto the accretion disk is expected to be lower for the pencil beam geometry than for the fan beam case and this could explain the change on the decay trend in the optical light curves from MJD ∼58112 onwards (see Fig. 6).Indeed, the secondary maxima in the pulse profiles starts shifting on phase from MJD 58112 to 58139 pulse profiles (see Fig. 1 of Liu et al. (2022a)).Liu et al. (2022a) observed that the pulse profile changed to double-peaked again at around MJD 58139 (dashed-dotted gray line in Fig. 6), and this change has been interpreted as the transition from super-to sub-critical accretion and as the offset of the accretion column and the presence of a pencil beam dominated accretion (Doroshenko et al. 2020;Liu et al. 2022a).This is in agreement with the fact from MJD 58140 onwards, the optical/UV SEDs can be reproduced just by the heating of the companion star, and the irradiation of an accretion disk is nos necessary to fit the observations.Another possible explanation to this optical hump be that a change in the geometry of the Be circumstellar disk had occurred.A change in the inclination and/or in the precession of the Be disk could be responsible for an increase of the optical/UV emission (with a major contribution at redder wavelengths), and could also be related with the end of the outburst.

Summary and conclusions
We have studied the optical/UV emission during the 2017-2018 Super-Eddington giant outburst of Swift J0243.6+6124.We have considered the possible origins of this emission and concluded that the main mechanism that can explain the observed light curves is thermal reprocessing of the X-ray emission by the donor star and the accretion disc.
We have estimated an interstellar reddening of E(B − V)=1.20±0.01mag from fitting the SED corresponding to X-ray quiescence.We have calculated the flux conversion factors to be applied when converting observed magnitudes in the UVOT filters to monochromatic fluxes.We have joined these UVOT fluxes with the compilation of all the V-Johnson band observations from Reig et al. (2020) to build the SED for each day of the outburst (whenever there were available observations).
We have studied the decay timescales of the optical and UV light curves which are between 2-4 times the X-ray light curve timescale.We have also calculated the correlation indexes of F opt/UV ∝ F β X and obtained values of β which depends on the considered X-ray band and on the optical/UV filter, and are larger for the rise than for the decay phase of the outburst for the redder wavelengths.This difference between the rise and decay is not observed for the UVOT-um2 and UVOT-uw2 filters.The correlation indexes are close to β ∼ 0.5 when considering bolometric Article number, page 14 of 16 X-ray luminosity during the rise of the outburst, which is compatible with the classical value expected from reprocessing proposed by van Paradijs & McClintock (1994), but is smaller when the soft X-ray energy bands are used and for the decay of the outburst.Then, to be able to reproduce the observed optical/UV fluxes of Swift J0243.6+6124, a more sophisticated approach is required.
For this reason we built a physical model including the contribution of the X-ray heating of the companion star surface, the irradiation of the accretion disk, and the contribution from a viscously heated accretion disk.As an input of the model, we estimated the X-ray luminosity from the spectral fits performed to all the available Swift/XRT data considering a power law plus a black body model, and then we added a high-E cutoff with variable folding energy to fit the Swift/BAT observed count rates.From the SED fits, we obtained that the principal mechanism contributing to the UV fluxes along the whole outburst is the heating of the companion star.We also concluded that the irradiation of the accretion disk needs to be considered in order to reproduce the observed optical fluxes SEDs from MJD 58047 to 58120 and contributes significantly to the UV emission close to the peak of the outburst.Before and after these times, the heating of the companion star is enough to reproduce the observed optical/UV fluxes, although the SEDs after MJD 58120 only have one point per day (mostly in the UV filters).
According to our modelling, the emission from the viscously heated disk does not contribute significantly to the optical and UV fluxes during the giant outburst of Swift J0243.6+6124.As a first attempt, we performed the SED fits, leaving the outer radius of the accretion disk R out as a free parameter.From the values providing the best fits, we observed that the fitted value of R out first increases until reaching value of R out ∼0.1 a which remains constant from ∼MJD 58047 until the peak of the outburst.After that, the fitted value of R out increases from MJD 58062 to MJD 58087 reaching R out ∼0.2 a, which is roughly the Roche Lobe size of the NS.After that, the fitted value of R out remains constant at this value.This could be interpreted as an increase of the size of the external radius of the accretion disk and this could happens either if the accretion disk grows or dissipates.A change on the factors considering changes in the geometry of the inner flow, the fraction of absorbed emission, and the possible anisotropy of irradiation could also reproduce a similar effect.Indeed, when fixing R out ∼0.2 a and leaving a combined factor including these effects (k) as a free parameter in the SED fits, the results are compatible with some of the scenarios proposed in previous works studying X-ray spectral and pulse profiles variability.
An alternative explanation to the apparent increase of the outer radius of the accretion disk or to the change in the inner flow accretion geometry after the peak of the outburst, would be that a change in the geometry of the circumstellar disk took place.Such a change could introduce an extra contribution (more significant at redder wavelengths), but our few observational data points do not allow to disentangle between both scenarios.
The mechanisms considered in the model presented in this work are suitable to explain qualitatively the observed slopes of the optical/UV SEDs and the variability along the outburst.Under the considered assumptions, they are able to reproduce the observed fluxes.In the case of a Be disk not co-planar to the orbit, the irradiation of the Be disk should play a role.In that hypothetical case, not including this mechanism could lead to an overestimation of the fitted values of R out , and would require to redefine some of the considered assumptions, such as the values of the albedo or thermalization efficiency, which could be dif-ferent than the considered ones and even variable along the outburst.Correct modelling of the Be disk structure and geometry would require complex theoretical considerations and additional observations, and such a detailed analysis is proposed for a future work.

Fig. 3 .Fig. 2 .Fig. 3 .Fig. 4 .
Fig.2.The normalized transmitted counts through the Swift/UVOT filters for Vega (green) and for the theoretical model of an O9.5 V star reddened with E(B-V)=1.20 mag (blue).They have been obtained multiplying the flux density by the effective area of each filter and converting fluxes to count rates.The effective wavelengths for each filter are marked with gray dashed lines.The largest differences are found for the UVOT-uw1 and UVOT-uw2 filters.

Fig. 5 .
Fig. 5. Fit of the long-term variations of the observed fluxes in the optical V-Johnson band.These estimated fluxes (in blue) have been removed from the observed fluxes (in red) to estimate the extra emission due to the outburst (see Sect. 3.3).

Fig. 6 .
Fig. 6.Optical vs X-ray bolometric (15-50 keV) emitted luminosities during the bright outburst of Swift J0243.6+6124.The optical V-band integrated fluxes have been obtained multiplying the monochromatic fluxes (converting magnitudes to fluxes with a zeropoint ZP=3.631×10 −9 erg s −1 cm − 2 Å −1 ) by the effective width of the filter (890Å) and dereddened considering E(B-V)=1.20 mag A distance d=5.2 kpc has been assumed in the luminosity calculations (L = 4πd 2 × F).The times related to changes in the pulse profiles and the X-ray spectrum proposed by Liu et al. (2022a) and Doroshenko et al. (2020) are marked with blue and pink lines respectively.The gray dashed-dotted line was proposed as a critical time (MJD 58139) on both works (see Sect. 4.3 for details).

Filter
Fig.7.Optical/UV vs X-ray (0.1-500 keV) logarithmic fluxes along the outburst of Swift J0243.6+6124.The correlations show different slopes during the outburst rise than during the decay epochs for the filters with longer wavelengths (V-Johnson to UVOT-uw1).In this case, different fits have been performed for the rise (blue line) than for the decay (red line).For the UVOT-um2 and UVOT-uw2 filters, a global fit (purple line) has been performed.

Fig. 8 .
Fig.8.Models considered in this work to fit the optical/UV observations (red points) and XRT spectra (purple line) for days 58062 MJD (up) and 58094 MJD (bottom).The contribution of X-ray heating of the companion star is represented with a pink dotted line, the irradiated disk is plotted as a brown dashed line, the viscously heated disk is plotted with a green dash-dotted line.The xspec models bbodyrad and pegpowerlaw used to fit the XRT spectra are plotted with light blue dotted and navy blue dashed lines respectively.

Fig. 9 .
Fig. 9. Illustration of the effect of X-ray heating of the surface of the companion star.The calculations for the peak of the outburst (58064 MJD) are shown.

58120Fig. 10 .
Fig. 10.Optical/UV SEDs evolution along the outburst.The reddening corrected observed fluxes are plotted with different colours: V-Johnson in red, UVOT-b in yellow, UVOT-u in green, UVOT-uw1 in blue, UVOT-um2 in navy blue, and UVOT-uw2 in purple.Error bars are plotted in gray.Best fitted model predictions are plotted for comparison (black solid lines).The pink dotted lines represent the modelled contribution by the extra emission from heating the surface of the companion star.The brown dashed lines correspond to the contribution from an irradiated accretion disk with the value of R out providing the best fit.Article number, page 13 of 16

Fig. 11 .
Fig. 11.Outer radius of the accretion disk which provided the best fit for each SED.

Fig. 12 .
Fig. 12. Results of the SED fit when the outer radius of the accretion disk is fixed to R out =0.2a, and we use as free parameter the factor correlated with the geometry and fraction of the radiation reprocessed (k).
Optical vs X-ray long-term light curves of Swift J0243.6+6124.the light curve contains observations from different sites: green points come from ASAS-SN Variable Stars Database, pink points correspond to converted Swift/UVOT-v observations, purple points are from AAVSO, orange points are from Skinakas observatory, and light blue points are from Aras de los Olmos Observatory (see Sect. 2.1 for details).The navy blue line represents the X-ray count rate light curve measured with Swift/BAT in the 15-50 keV band.