Performance and first measurements of the MAGIC Stellar Intensity Interferometer

In recent years, a new generation of optical intensity interferometers has emerged, leveraging the existing infrastructure of Imaging Atmospheric Cherenkov Telescopes (IACTs). The MAGIC telescopes host the MAGIC-SII system (Stellar Intensity Interferometer), implemented to investigate the feasibility and potential of this technique on IACTs. After the first successful measurements in 2019, the system was upgraded and now features a real-time, dead-time-free, 4-channel, GPU-based correlator. These hardware modifications allow seamless transitions between MAGIC's standard very-high-energy gamma-ray observations and optical interferometry measurements within seconds. We establish the feasibility and potential of employing IACTs as competitive optical Intensity Interferometers with minimal hardware adjustments. The measurement of a total of 22 stellar diameters are reported, 9 corresponding to reference stars with previous comparable measurements, and 13 with no prior measurements. A prospective implementation involving telescopes from the forthcoming Cherenkov Telescope Array Observatory's northern hemisphere array, such as the first prototype of its Large-Sized Telescopes, LST-1, is technically viable. This integration would significantly enhance the sensitivity of the current system and broaden the UV-plane coverage. This advancement would enable the system to achieve competitive sensitivity with the current generation of long-baseline optical interferometers over blue wavelengths.


INTRODUCTION
The interferometry technique combines visible light (or other electromagnetic wavelengths) employing two or more telescopes to obtain a higher resolution image by applying the principle of superposition.Radio interferometers in effect record the electromagnetic field at disparate locations, and from these data an aperture, which can be as large as the Earth (e.g., Event Horizon Telescope Collaboration: Akiyama et al. 2022), is synthesised via software.In the optical wavelength domain, achieving similar resolution with kilometerscale aperture synthesis is conceivable.However, direct recording of the electromagnetic field in the same manner as in radio is currently not feasible in optical interferometry.Consequently, optical interferometry resorts to employing more indirect methods.
The best-known technique for optical interferometry goes back to Michelson & Pease (1921) and involves bringing the light from different telescopes together and making it interfere.The challenge is to maintain coherent optical paths between the telescopes.The early work could do so over separations of about 10 m, and longer coherent baselines would not be realized for several decades.
In between, however, Hanbury Brown & Twiss (1958) developed a different kind of optical interferometry technique named intensity interferometry, which involves correlating intensity fluctuations at different telescopes, without the need of physically combining optical beams, and not requiring coherent baselines. 1The Narrabri Stellar Intensity Interferometer (NSII) extended the then-new Hanbury Brown and Twiss (HBT) technique to baselines between 10 and 188 m.Over the 1960s and early 1970s the NSII observed 32 individual stars and 9 multiple-star systems, measuring stellar diameters and astrometry of binary systems (Hanbury Brown et al. 1974b).Its limitation was not the baseline but the signal-to-noise ratio (S/N), which prevented observing fainter sources.
Interest then returned to the Michelson-Pease type of optical interferometers, which eventually also achieved baselines of hundreds of metres, without the S/N limitations of HBT interferometry.CHARA and VLTI are the best known of these, and have had many successes, which we will briefly discuss below.Extensions to km-scale baselines are under development (see e.g., Bourdarot et al. 2020) with even more ambitious proposals involving going to space or the Moon (Labeyrie 2021), but none of these projects are imminent.
For optical interferometry at km-scale baselines a nearer-term prospect is intensity interferometry with new generation instrumentation.In particular, implementation as a second observing mode in Imaging Atmospheric Cherenkov Telescopes (IACTs) has been advocated for some time (Le Bohec & Holder 2006;Dravins et al. 2013), and the last few years have seen a revival of intensity interferometry.Early results from MAGIC have been reported in Acciari & al. (2020b) and Cortina et al. (2022).VERITAS, following the initial implementation of intensity interferometry (Abeysekara & al. 2020), has announced an ongoing survey of stellar angular diameters (Kieda et al. 2022).H.E.S.S. has developed intensity interferometer (Karl et al. 2022).In addition to IACTs, results from intensity interferometry through standard optical telescopes have also been reported (e.g., Horch et al. 2022;Matthews et al. 2023).
The recent scientific results from HBT interferometry have been mainly stellar radii and their implications.These are modest compared to what Michelson-Pease interferometry has achieved over recent years (Eisenhauer et al. 2023).Here we list some examples of recent scientific achievements led by the technique, with a special focus on those topics in which intensity interferometers may contribute: 1) Radial oscillations of Cepheids have been spatially resolved by CHARA (e.g., Nardetto et al. 2016) and intensity interfometry has the potential to resolve these and more complicated asteroseismic modes.2) For fast-rotating stars, CHARA has resolved (e.g., Che et al. 2011) rotational flattening accompanied by gravity darkening of the equator compared to the poles.As shown by Nuñez & Domiciano de Souza (2015), the simultaneously available UV coverage provided by intensity interferometry telescope arrays is ideal for this science case.3) The surface of Betelgeuse has been resolved with the VLTI (Montargès et al. 2021).Image reconstruction algorithms have also been adapted to intensity interferometry (Dravins et al. 2012), and a combination of both techniques would lead to a wider resolution coverage.4) Especially well known is the astrometry of relativistically moving stars in the Galactic-centre region using VLTI (e.g., GRAVITY Collaboration: Abuter et al. 2022).5) More ambitiously, gravitational-wave-emitting binaries have also been discussed as possible targets for intensity interferometry (Baumgartner et al. 2020).
In order to demonstrate the viability of these scientific objectives and justify farther pursue of this technique, the performance of the current instrumentation needs to be studied in more detail, as well as exploring the systematics that may affect instruments far exceeding the sensitivity of the NSII.

HARDWARE SETUP
MAGIC is a system of two IACTs located at the Roque de los Muchachos Observatory on the island of La Palma in Spain (Aleksić et al. 2016a).Equipped with 17 m diameter parabolic reflectors and fast photomultiplier (PMT) Cherenkov-imaging cameras, the telescopes record images of extensive air showers in stereoscopic mode, enabling the observation of very-high-energy (VHE) gammaray sources at energies of few tens of GeV up to tens of TeV (Aleksić et al. 2016b).
In April 2019, to prove that MAGIC was technically ready to perform intensity interferometry observations, a test was performed using MAGIC telescopes and an oscilloscope as a readout (Acciari & al. 2020b).Temporal correlation was detected for three different stars of known angular diameter.The sensitivity and the degree of correlation were consistent with the stellar diameters and the expected instrumental parameters.However, the acquisition was constrained by a low duty cycle, and the mechanical installation of necessary optical filters hindered a swift transition between VHE and interferometry observations.
In the following we describe the technical modifications that have been implemented in MAGIC to enable intensity interferometry observations, summarised with a diagram in Fig. 1.Some of these modifications have already been discussed in detail in Acciari & al. (2020a); Delgado & al. (2021); Cortina et al. (2022).A key consideration was to not affect regular VHE observations, allowing a smooth and effortless transition between "VHE observation mode" to "interferometry observation mode" and back in less than one minute.Since the implementation of these modifications in 2021, a total of 192 observing hours have been performed between January and December of 2022, predominantly during bright Moon-light periods.

Photon detectors and signal transmitters
MAGIC telescopes are equipped with 1039-pixel PMT cameras at their primary focus.The PMTs are 25.4 mm in diameter and have 6 dynodes.A hexagonal shape Winston Cone is mounted on top on each PMT.The distance between PMT centers is 30 mm and corresponds to a 0.1 • FOV.
The camera and data acquisition of the MAGIC telescopes were developed for efficient detection, recording, and offline reconstruction of temporally brief VHE gamma-ray events (Bitossi et al. 2007).The detection of low-energy events on a ns time scale against the constant night sky background (NSB) requires a wideband, low dispersion signal path from the camera focal plane to the data acquisition trigger and readout.The camera uses fast PMTs with high quantum efficiency (QE) in the wavelength range of incident Cherenkov light to respond to events while being less sensitive at NSB wavelengths.A dedicated hardware trigger allows the extended 2D image data to selectively record only when signal is present (Dazzi et al. 2021).A high dynamic range of 10 3 supports energy reconstruction over a wide range of incident energies.Slow control monitoring of several operating parameters including the PMT anode direct currents (DC) is done continuously during observations.As a summary, the specifications of the camera and data acquisition for VHE observations are: • High detector QE at Cherenkov wavelengths, suppressing NSB.
• Wideband, low dispersion signal path from the camera to the separate readout.
• Low noise in the signal path to facilitate trigger response to few-photon signals.
• Isochronicity at the ns level across all channels.
• Selective recording of pixelized 2D Cherenkov images with a dynamic range of ∼ 10 3 .
• Continuous PMT DC monitoring at one second intervals (for monitoring purposes).
In contrast, SII observations require a continuous recording of optical signals from a point-like source over observation times ranging from minutes to hours, with a different set of desirable characteristics: • High detector QE in the wavelength band of interest.
• Wideband, low dispersion signal path from the camera to the separate readout.
• Low noise in the signal path compared to the total starlight signal.
• Very low correlated noise in the signal path.
• Continuous recording of a few pixels per telescope with moderate dynamic range and fixed temporal delay between channels (low jitter).
• Continuous PMT DC monitoring at one second or less intervals.
The dynamic range necessary to sufficiently resolve the steadystate input signal is dictated by the expected photon flux divided by the duration of each photon signal.This dynamic range is typically more than an order of magnitude less than that implemented for VHE observations.To aggregate measurements across varying flux levels, an essential prerequisite is a quantification of the input photon flux.In this context, we rely on the slow control DC reports, which encompass system logs storing information on PMT average photon currents and the applied high voltage (HV) over time.
Aggregation of measurements at different flux levels requires a measure of the input photon flux, for which we use the slow control DC reports (available system logs storing PMT average photon currents and applied high voltage (HV) as a function of time).
Instead of having additional photo-detectors mounted on top of their cameras (as done by VERITAS or H.E.S.S.), MAGIC-SII observations are performed with the same pixels employed for gammaray observations.Only a few selected camera pixels and their optical analog signal transmission are used in conjunction with a separate signal receiver and digitizer optimized for SII in place of the VHE readout.The used camera pixels contain a PMT and wideband AC coupled amplifier which drives vertical cavity surface emitting lasers (VCSELs) operating at 850 nm for analog optical signal transmission.The characteristics of the camera pixels and their optical transmission are detailed in Aleksić et al. (2012).
The SII receiver consists of a multimode fiber coupled, battery reverse-biased photodiode with responsivity suitable for 850 nm and a bandwidth of 4 GHz.The detector is placed in an radio frequency (RF) enclosure coupled to the input of a Femto HSA-Y-2-40 wideband amplifier which supplies an AC coupled output signal level of ∼ 15 mV amplitude per single phe (dependent upon PMT gain settings)  to a 50 Ω semi-rigid coaxial connection to the digitizer.The Femto amplifier has two separate output buffers, facilitating simultaneous operation of two different digitizers or concurrent signal monitoring.
The noise contributions originating from the electronics can be classified in two different categories depending on their contribution to the resulting correlation and its uncertainties: • Uncorrelated noise: The VCSEL used in analog optical transmission contributes as relative intensity noise (RIN).Since the signal variation is small compared to the bias current of the VCSEL, it can be considered to be fixed and independent of the input signal level.The VCSEL RIN is also subject to random variations in amplitude and spectrum.The pixel preamplifier and Femto amplifier used in the receiver both contribute with additional wideband uncorrelated noise to the output signal.This noise contribution is of fixed amplitude and is smaller than the VCSEL contribution.Their combined contribution can be disregarded in the analysis provided the total optical signal power per measurement interval is greater than the total uncorrelated noise power from the electronics, otherwise it needs to be dealt with as a correction term when calculating the uncertainty.Other sources of uncorrelated noise may be pickup local to a particular telescope camera or an individual digitizer channel.
• Correlated noise: Channel to channel crosstalk is an example of delay-dependent correlated noise, which is heavily correlated near delay zero between channels and uncorrelated otherwise.In our setup we avoid crosstalk contributions by introducing a fixed optical delay on the order of 500 ns in adjacent digitization channels to locate the signal outside the region of disturbance.The SII optical receivers use RF shielding, battery photodiode bias, separate amplifier power adapters, and coaxial connecting cables with effective shielding to suppress crosstalk between receivers.External signals due to a source common to multiple channels such as Radio Frequency Interference (RFI) pickup in a common camera or adjacent computing equipment, may be transient or continuous in nature and independent of the delay between channels.Their contribution may render the anal-ysis ambiguous whenever they become comparable to the analysis uncertainty or signal amplitude.
The pixel preamplifier, optical signal transmission, and optical receiver have a combined transfer function magnitude bandwidth measured to be > 400 MHz with low dispersion using an R+S ZVA-8 vector network analyzer.Due to this, the signal channel bandwidth will be mostly determined by the PMT pulse time response (∼ 2 ns FWHM) and the antialiasing filter of the digitizer (125 MHz).The aggregate bandwidth of this setup is fairly consistent with the one used by Acciari & al. (2020b), approximately 110 MHz.Improving the bandwidth of the system will be discussed in section 4.5.
The DC measurement branch is used as implemented for monitoring purposes in the MAGIC cameras.Since its properties are critical to an accurate determination of the relative incident photon flux, a detailed description is provided.As shown in Fig. 2, the DC monitoring path consists of a resistive divider connected directly between the PMT anode and signal ground potential.The divider decouples the DC monitoring path from excessive PMT output levels and isolates the RF signal path from the monitoring circuit.The PMT current is sensed as a (negative) voltage at the output of the divider, which is amplified by a low offset, low temperature and time drift, low input bias current operational amplifier (OPA335, Burr -Brown from Texas Instruments) in inverting configuration.
The output of the sensing circuit is connected to a 12-bit ADC with internal reference and temperature compensation, MAX1231 from Analog Devices (formerly MAXIM IC), and read out via serial interface at 1 second intervals.The reference ground for all monitored values is separate from the pixel power return line to preclude voltage drops along the connection lines.

Optical filters mounts
As described by Hanbury Brown et al. (1974b), the S/N of the correlation of telescope signals is insensitive to the width of the optical Exclusively the effect of the angle distribution is shown here.Reflectivity and optical performance are expected to modify the transmission, but not the shape of the optical passband.passband of the detected light.However, a filter spectral response with sharp spectral cutoffs does improve the sensitivity of the measurement.In addition, accounting for MAGIC photo-detection efficiency and gain, we are observing bright stars capable of damaging PMTs after a short exposure time, given the excessive current flowing through the dynode system.For these reasons interferometry observations require installing narrow-band optical filters in front of the PMTs connected to the correlator.We are using interference filters manufactured by Semrock of model 425-26 nm.The spectral transmission curve for incident parallel light is centered at 425 nm and has a FWHM of 26 nm with relatively sharp edges.This shape is strongly modified in our setup because MAGIC has a low f/D ratio (close to 1).We have used the MyLight modeling online tool provided by Semrock2 to calculate the effective spectral transmission curve in a cone of half angle 26.5 • , as shown in Fig. 3.
The filters have a diameter of 50 mm and are held 20 mm in front of MAGIC photo-detectors using a mechanical frame ("filter holder").This frame piggybacks on an existing mechanical structure that holds a diffusely reflecting white plate used for absolute calibration of the reflectance of the mirror (the so-called "white target").This structure is operated remotely and can be deployed in front of the PMTs in a matter of seconds.Fig. 4 shows a picture of the filter holder under the white target.There is room for six filters in a horizontal line, all centered in front of PMTs that may eventually be equipped for interferometry.The additional filter slots also allow for a simultaneous signal and background monitoring using identical filters during interferometry observations.Each filter is placed and centered in front of the central PMT of the 7-pixel clusters of the MAGIC camera.The pixels used in this work have a ∼ 24 cm offset (0.8 • on sky) with respect to the center of the camera.

Active mirror control
The parabolic-shape reflector of a MAGIC telescope is approximated by using 246 individual mirror facets of spherical shape.The facets are fixed on the reflector according to their radius of curvature, which varies between 34 to 36 m.Eleven groups of such mirrors lead to a very low time spread of synchronous light pulses (< 1 ns) (Bastieri et al. 2008;Aleksić et al. 2012).The light-weight design of the MAGIC reflector, made of reinforced carbon fiber tubes, allows for fast slewing of the telescopes but is subject to deformations of the dish and sagging of the camera.Each mirror facet of 1 m 2 is equipped with 2 actuators, correcting for these effects via an Active Mirror Control (AMC) system (Biland et al. 2008).
The flexibility given by the AMC has been a key asset to the MAGIC intensity interferometry setup.As described by Hanbury Brown et al. (1974b), the sensitivity of SII strongly relies on the photon detection capabilities of the telescopes.The MAGIC Collaboration devotes significant observational, scheduling, and analysis efforts for keeping the bending models of the telescopes up-to-date to ensure that the AMC minimizes the effect of structural deformations of the telescopes, maximizing the amount of starlight reaching the right camera pixels (Wallace 1994).During interferometry-mode observations, the on-site crew performs a re-focusing of the AMC whenenever the zenith distance changes by at least 5 deg, to ensure the right bending model is applied as stars drift during the observation.
In addition to correcting deformations, the AMC allows for a broad range of configurations, transforming the MAGIC telescopes into a very versatile optical interferometer: • Full-mirror: as described in section 2.2, light from the target star focuses into a single pixel located with a 0.8 deg offset with respect to the camera center.This offset adds negligible spread to synchronous light signals (still < 1 ns).This is the standard observing mode of the MAGIC-SII setup, which focuses starlight into any of the 6 pixels behind the optical filter holders.
• Chess-board: by focusing half of the mirrors to one interferometry pixel and the other half to another (in a pattern similar to a chess-board), each MAGIC telescope becomes 2 virtual overlapping telescopes.As described in section 2.4, the GPU correlator is able to handle 4 input signals (computing 6 correlations), so we are able to simultaneously gather long-baseline correlation signals as well as a measurement close to our zero-baseline correlation.
• Sub-reflectors: as any combination of the mirrors located within MAGIC reflectors may be focused to the interferometry pixels, the MAGIC-SII system has the capability of sampling short-baseline correlations within the 1-17 m range by creating multiple sub-reflectors (of e.g.3-5 m in diameter).Even if the number of feasible targets to be studied with this configuration is limited (very-bright, very-large stars), the flexibility on the number of sub-reflectors and their relative location allows a broad sampling of the Fourier space of the image (UV coverage).This could be expanded to more complex setups, such as the I3T concept, envisioned by Gori et al. (2021).

Digitizer and correlator
The correlator hardware and software have been designed to harness the massively parallel nature of state-of-the-art GPUs to process in real time the data captured using two fast digitizer boards Spectrum M4i.4450-x8 PCIe 2.0.
Each digitizer deals with two channels providing up to 500 MS/s of simultaneous sampling rate with a resolution of 14 bits per sample.Furthermore, the selected digitizers support Remote Direct Memory Access (RDMA), enabling direct data transfers between the digitizer's memory and the GPU's memory.This eliminates the need for intermediate copies, resulting in increased throughput and reduced latency.The two Spectrum boards clocks are synchronized by means of the Star-hub module attached to the carrier card in the correlator chassis.Because using the High Frequency mode introduces a strong flip-flop correlated noise into the data stream, buffered mode (with a resistance of 50 Ω) is used for the digitizer (imposing the ∼130 MHz bandwidth described in the previous section).
The correlator is implemented in a computing server with off-theshelf hardware: two processors (20 cores in total), Solid State Drives (SSDs) for fast access and Hard Drive Disks (HDDs) for longer term storage and tests, and a Nvidia Tesla V100 GPU.The GPU chosen is the PCIe 3.0 x16 model with 5120 cores, 32 GB HBM2 RAM and 14 TFLOPs of single-precision performance.
The correlation code, described in section 3.1, has been used to process data from the two digitizers in real time at 4 GB/s.Alternatively, the code can dump raw data from one digitizer (two channels) to disk in real time for short periods of time, mainly for testing purposes.

DATA ANALYSIS
We use the prescription described by Acciari & al. (2020b) to perform intensity interferometry astrophysical measurements with the MAGIC-SII setup, where we define the contrast , proportional to the squared visibility  2 : where  is a constant,  is the Pearson's correlation at the  0 where the signal is expected,   are factors applied to correct for changes in the HV of pixel i, and   are the DC measured in pixel i, which is proportional to the photon flux.This expression assumes the signal is dominated by the stellar flux, and ignores the contribution of the NSB.When NSB currents become a significant fraction of the measured signal, this expression needs to be expanded to account for the contribution of the un-correlated photons: where  accounts for the ratio between stellar and NSB photon fluxes: Squared visibility measurements are calculated using eq.( 2) by performing the following steps: (i) () computation by the correlator, (ii) calibrate () with   ,   and  factors, (iii) variable time delay  correction, (iv) extract  at the expected  0 .
In this section we will describe the implementation developed for the low-level analysis, i.e. the GPU-based online computation of Pearson's correlation (), signal calibration and high-level analysis of square visibility measurements.

GPU-based on-line correlation
The main functionality of the correlator software is to compute the Pearson's correlation (()) for all possible pairs from the 4 channels used as a function of the time shift  between them over a window wide enough to cover the expected range of delays (in our case, a generous ± 2048 ns).It also computes the autocorrelation of these channels, resulting in a total of 6 correlations and 4 autocorrelations.To perform these computations efficiently, the software exploits the convolution theorem in frequency domain implemented with Fast-Fourier-Transform (FFT).The correlation for each pair in frequency domain is computed for time frames of a given size (generally 2 18 samples per channel), added in sets of 500 frames and then converted back into  domain using the inverse FFT.The resulting correlations, normalized with statistics of the input signal to obtain (), are time stamped and stored in disk together with the mean voltage and standard deviation of each input channel used for the computation.
The software is developed in CUDA C using the Nvidia FFT library3 for the convolution computation, and the Spectrum CUDA Access for Parallel Processing (SCAPP) SDK 4 for transferring the data between the digitizers and the GPU in streaming mode.It is divided in three parts: (i) The initialisation section configures several parameters of the digitizers and the GPU (e.g., number of channels, number of correlations, acquisition rates, input paths and ranges, RDMA, execution time), initialises the data structure in both the CPU and the GPU and creates the thread that writes the resulting correlation to a storage media.The data processing loop then starts using a double buffering scheme for the input and output, which allows for a live-time of the correlator of ∼100%.
(ii) The data processing loop runs continuously for the programmed duration of data taking (generally, a data run of 5 minutes), until an error occurs or until it is interrupted by the operator.In this loop: (a) The code first checks and retrieves into an input buffer a time frame of a given number of samples from each stream of incoming data from the channels of the digitizers.
(b) Statistics, such as the mean and the standard deviation, are obtained from the input buffer and added to accumulators in an output buffer for each input channel.
(c) Simultaneously, the FFT is computed for the data in the input buffer of each channel.
(d) After this, cross correlations and auto correlations are computed by multiplying the obtained FFT and complex conjugate for each channel, which are then added to an accumulator in the output buffer.
(iii) Every 500 iterations in the data processing loop, the output buffer is swapped and a writer thread is awoken.This thread normalizes (using the computed statistics accordingly to the definition of the Pearson's correlation) and saves the accumulated correlation results and any obtained statistics to disk in binary format, sets the accumulators to zero and goes back to sleep.
In addition to its main functionality, the software can be used to store the raw signal of one digitizer directly to disk.In this mode a sample with a resolution of 14bits is saved every 2 ns for each channel.Due to the size of raw data (1 Gb/s/channel), only data from two channels can be saved simultaneously, typically in runs of 10 s.

Coherence calibration and signal extraction
The correlator described in the previous section computes and stores () in a wide range of time delays (± 2048 ns).One () array is calculated from a fixed integration time, generally lasting ∼ 250 ms, and stored with their corresponding time tag.As introduced in this section, these signals need to be calibrated using the   ,   and  factors, as shown in eq. ( 2).As the time  0 in which the correlation signal is expected changes along the observation (as the star moves across the sky), () arrays need to be aligned in time with respect to the expected  0 before averaging.
The MAGIC SII setup employs a total of 3 pixels in each telescope.As shown in Fig. 4, multiple filter holders are available.Currently, a total of 3 Semrock 425-26 nm optical filters are installed in each telescope: two pixels in each telescope, tabulated in Table 1, are used for signal extraction (connected to the correlator) while a third pixel is used exclusively for simultaneous background DC measurements over the same optical passband.DC reports from the telescope slow control system store the DC and HV of each pixel once per second for each camera.DC conversion factors are calculated from calibration measurements, allowing to calculate ( )  at any time, and therefore also ()  and  terms.Most interferometry observations are performed with the same HV values, but those that were performed with a different configuration are still usable as long as their gains are corrected with the   terms, also inferred from calibration observations.The time delay correction is applied to each () array stored by the correlator, by subtracting the expected time-delay of the correlation signal  0 , so that the expected correlation signal will always be located at  ∼ 0. The expected time-delay  0 depends on two terms: 1) the specific hardware delay from each channel used in the observation (dominated by the length of the optical fibers) and 2) the time of flight difference between photons detected in each telescope.As current digitizers operate at 500 MHz, arrays are aligned in steps of 2 ns.
Once the correlation signal of a given source has been calibrated and time-delay corrected (i.e.

𝜌(𝜏− 𝜏
) it is collected and accumulated, generally in bins of baseline.The averaging is performed by weighting each array by the measured off-peak variance (using values far from  ∼ 0).As the uncertainty expected of the correlation is inversely proportional to the photon flux measured, a linear average across time would magnify the uncertainty in the correlation when combining observations with different quality (measuring different fluxes, due to cloudiness, pointing direction or hardware issues).After this weighted average, contributions of bad-quality data are minimized according to their expected larger variance, allowing to combine observations with very different data quality with minimal negative impact on signal to noise.
After the correlation is averaged in time, the only step necessary to compute , as shown in eq. ( 2), is to extract the amplitude of the correlation signal.This is done by fitting a Gaussian function, where the resulting amplitude will be  (see Fig. 5).Given the limited bandwidth of the digitizer, data is strongly correlated, and therefore its bi-dimensional correlation matrix is used in the  2 minimization to properly weight the residuals.As will be justified in section 4.4, a high-pass 12 MHz cut filter is applied (both to the data and the Gaussian fitted function) to mitigate a faint correlated noise seen when integrating over long observing periods.The uncertainty of the squared visibility measurement Δ is set to the standard deviation of the off-peak correlation in the time-delay range surrounding the expected correlation signal, i.e., from -140 to -40 and from 40 to 140 ns in  −  0 .
MAGIC PMT gains are calibrated during VHE gamma-ray observations by using interleaved calibration pulses Aleksić et al. (2016a,b).By digitizing with the GPU data acquisition (DAQ) these same calibration pulses we are able to store the average response of these pixels to bright short light pulses.Table 1 shows the width of calibration pulses measured by the different pixels equipped for interferometry observations.In addition, as the GPU DAQ also computes the auto-correlation of each channel, we can ensure starlight observations are providing the expected bandwidth.From the average response to calibration pulses between different pixels we are able to calculate the expected shape of our correlation signal ( ∼ 2.2 ns).This value already includes the fixed jitter (0.4 ns) we expect for full-mirror observations from MAGIC I (the first telescope built in 2004) since mirrors in the reflecting dish were installed in two layers of facets separated by 60 mm in the optical axis direction.As deviations with respect to a Gaussian are at a few percent level and statistical uncertainties from individual correlation signals will not reach such precision, we use a Gaussian function to extract the amplitude of the correlation signal.The width of the measured pulses is currently limited by the bandwidth of our digitizers., from ∼ 30 min of observation with an average baseline of 52 m pointing to eps CMa (Adhara).The amplitude of the correlation is shown in the legend, using a Gaussian function with a 12 MHz cut high-pass filter applied.Best fit is drawn in blue while 1- uncertainty region is drawn in orange.Bottom) calibrated DC (in A) from each MAGIC telescope as a function of the zenith distance of the observation, both for the current associated to the star, as well as the NSB.Different DC values from the same telescope and zenith distance come from observations from different nights.Wind gusts produce relatively large fluctuations in M2 (less protected from wind than M1) while the low-amplitude sawtooth effect seen in both telescopes is due to the azimuth tracking.

Angular diameter analysis
As described in section 3.2, correlation data is calibrated, timedelay corrected and bundled to perform the weighted average over the observed time.In this work, data is bundled in uniform steps in baseline, and  2 measurements are extracted from the amplitude of each correlation signal.In the following we will consider either uniform disc (UD) or limb darkened (LD) profiles for the observed stars, i.e. radially symmetric models, therefore baseline  is the only relevant parameter to consider.In the UD scenario, visibility  can be expressed with the Bessel function of the first order  1 , as in e.g.Berger & Segransan (2007): where   is the diameter of the uniform disc,  the projected inter-telescope distance (baseline) and  the central wavelength of the optical bandpass of our setup (420 nm, see Fig 3).When including the limb darkening effect, we follow the prescription introduced by Hanbury Brown et al. (1974a), described as ( 5) where   =   /,   is the limb darkening coefficient and   is the angular diameter of the limb-darkened star.
Stellar diameters are measured by fitting eq. ( 4) to the  2 () measurements.The statistical uncertainty of fitted parameters, such as Δ and the zero-baseline correlation Δ(0), are considered to be the largest value between the one extracted from the  2 minimization method (the value increasing the  2 by  2  , as in Newville et al. ( 2016)) or the one calculated via bootstrapping.As the size of the MAGIC reflectors is comparable to the distance between the telescopes, eq. ( 4) is evaluated following the true distribution of distances between random points within the two reflectors (see Fig. 6).
Under these assumptions, in order to constrain the diameter of stars with uncertainties down to the few percent level it is necessary to reach such signal to noise both in the  2 () and  2 (0).As discussed in Hanbury Brown (1974), the zero-baseline correlation of an intensity interferometer is a constant of the system that mainly depends on the electronic and optical bandwidth of the hardware setup.Given the MAGIC interferometer is only composed by a single pair of telescopes, only stars located in favourable locations in the rotating sky provide the possibility of measuring a wide range of baselines (in the 20 to 87 m range), to reasonably constrain not only , but also  2 (0).As the latter is expected to be constant for a given hardware setup, as soon as  2 (0) is properly constrained, its measurement will be used to measure  of stars in less favourable positions in the sky, which only allows for the measurement of a very limited range of baselines.For the MAGIC telescopes, the ideal star for constraining  2 (0) is eps CMa (i.e.Adhara).As shown in Fig. 5 and 6, eps CMa is bright enough to provide strong correlation signals over reasonable observing times, and allows the broad baseline coverage required to strongly constrain both  and  2 (0).
As described in Acciari & al. (2020b), because the signals from MAGIC are not DC coupled, other parameters may affect the measured  2 (0) value, such as gains of the PMTs, or those associated to the DC measurements.For this reason, from here on we assume independent zero-baseline correlation measurements for each pixel pair:  2 (0) −  .By performing a simultaneous  2 minimization to all available observations from eps CMa (with the 4 different channelpair combinations) and assuming a common , we improve the statistical uncertainty of those  2 (0) −  values with significantly less data (mainly A-D and B-C pairs).Note the UV coverage of all channel pairs is almost identical, so if the source had non-radially-symmetric features, the expected systematic added by this assumption is negligible.The stability of  2 (0) −  will be discussed in section 4.4, along with all the sources of systematics evaluated in this work.

RESULTS
In this section, we summarise some of the results we have achieved with on-sky observations with the MAGIC-SII system, covering its calibration, validation, and new astrophysical measurements of stellar diameters.Since the prototype hardware implementation used in Acciari & al. (2020b), the upgraded setup presented here was designed, installed, and commissioned by January 2022.Since then, in the period between January and December 2022, 192 hours of data have been acquired, using a diverse set of observing modes (see section 2.3 for a description of the interferometry observing modes): ∼ 101 hours of full-mirror observations using pixels 251 (labeled A-C, as in Table 1), 28 hours of chess-board observations (correlating the 4 combinations of A, B, C and D channels) and 63.4 hours of full-mirror observations using pixels 260 (labeled B-D).
The targets observed over this period can be broadly classified into two categories: reference and candidate stars.We consider a target to be a reference star if their diameter has already been directly measured by other instruments over similar wavelengths (400-440 nm).We used the following selection criteria: angular diameter and declination allowing the determination of their stellar diameter with MAGIC baseline, as well as bright enough to detect correlation signals over reasonable observing times.From this selection of stars, we excluded fast rotators, spectroscopic binaries, those having bright close companions (Δ < 2.2 and distance below 10 mas) and those showing large variability in magnitude (Δ > 1).By measuring the diameter of reference stars we intend to confirm the validity of our analysis and hardware setup.Candidate stars are those not having a direct measurement of their diameter over similar wavelengths, but their predicted size and declination as well as their brightness allow a direct detection of their diameter.Their predicted diameter is extracted from several sources (Bourges et al. 2017;Swihart et al. 2017;Bonneau et al. 2006Bonneau et al. , 2011)).
All data was acquired and analyzed following the steps described in section 3.For simplicity, a fixed step size of 5 m in baseline was used for all sources, even if this binning size may not be optimal for the faintest stars presented.The significant amount of data presented in this work allow us to evaluate the performance of the system, to confirm the validity of analysis as well as evaluating the scale of the systematics associated to these observations.

Data and analysis consistency
Using the analyzed data we can test the scale of the residuals of  −  0 of the measured correlated signals.We select the correlated signals (each resulting from a set of observations bundled over constant steps in baseline of up to 5 meters) with a signal to noise larger than 3.No matter which pair of channels is used for the correlation, the location of all correlation signals fall within a ± 1 ns window (as shown in the left panel of Fig. 7, most within ± 0.5 ns from  −  0 = 0).The stability of the location of the correlation signal validates the locations of the telescopes assumed in the analysis, the hardware delays assumed for each channel as well as the calculation of photons time of flight as a function of pointing direction.
As discussed in section 2.1, the expected electronic bandwidth of the system is ∼ 125 MHz, dominated by the currently used digitizers.By fitting the width of the measured correlated signals we can test that the performance of the system matches expectations.As introduced in section 3.1, in addition to the 6 possible cross-correlations, the GPUbased correlator is able to store the auto-correlation of each channel.This provides an independent and simultaneous measurement of the electronic bandwidth of each channel (not including all possible effects, such as jitters between different pixels).As shown in the right panel of Fig. 7, the width of all correlated signals is consistent with a Gaussian sigma around 2.2 ns, consistent with the expectations computed both from calibration pulses (see section 3.2) as well as from measured autocorrelations.This excludes the presence of strong jitters between pixels, and shows that all pixels equipped for SII observations show very comparable performance.
The last consistency check performed to evaluate the validity of our analysis is to reconstruct all stellar diameters treating each channel pair as an independent dataset.As introduced in this chapter, the largest fraction of data employs A-C and B-D channels (101 and 63 h respectively) but an additional 28 h were taken in chessboard mode, in which all four correlations are possible (A-C, A-D, B-C, and B-D) although with half of the mirror area.In order to perform these channel-wise stellar diameter measurements,  2 ,  (0) are computed independently for each channel by performing a simultaneous fit of all sources with available  2 ,  .As shown in Fig. 8, stellar diameter measurements over the different correlation pairs show remarkable agreement, down to uncertainties in the few-percent level for at least two stars.

Stellar diameter measurements
Here we report stellar diameter measurements performed between January and December 2022.As introduced in section 4, only a few of the stars MAGIC-SII observed have a previous direct measurement of their angular diameter over similar blue wavelengths, allowing a one-to-one comparison of the measurements.The angular diameter measured by MAGIC as a function of the one measured by other instruments is shown in Fig. 9.The measurements performed by the MAGIC-SII system are consistent with previous measurements, both from the NSII, VERITAS and CHARA (Hanbury Brown et al. 1974b;Abeysekara & al. 2020;Gordon et al. 2019).All the stellar diameter measurements and physical parameters used to determine linear limb darkening coefficients (  ) together with the reference values we compare with are listed in Table 2.
In addition to the reference stars we observed, we selected 13 stars   Brown et al. 1974b), CHARA (Gordon et al. 2019;Challouf et al. 2014) and VERITAS (Abeysekara & al. 2020)).Physical parameters (effective temperatures and surface gravity) extracted from: Anderson & Francis (2012); Baines et al. (2018); Soubiran et al. (2016).Limb darkening linear coefficient   interpolated from Claret & Bloemen (2011) as done in Abeysekara & al. (2020).Measurements performed with the MAGIC-SII system are shown, assuming both uniform disc and limb darkened profiles, showing statistical uncertainties and maximum expected systematic deviation.for which their diameters have not been directly measured before in our bandwidth (400-440nm).Table 3 lists MAGIC new stellar diameter measurements (both using UD and LD models) together with the physical parameters used to determine   .They are mostly early-type stars between 2.07 and 3.73 magnitudes in B and estimated angular diameters between 0.3 and 1.3 mas.As shown in Fig. 10, MAGIC measurements are again nicely consistent with expectations.A significant fraction of these sources are known fast rotators, and due to the equatorial bulge produced by the centrifugal force they may deviate significantly from the radially symmetric models assumed here.In future works, the MAGIC Collaboration will release interferometric observations following community standards such as OIFITS data products (Duvert et al. 2017

Sensitivity evaluation
As introduced in Acciari & al. (2020b), from equation 5.17 in Hanbury Brown (1974) we are able to calculate the S/N we expect for a given correlation signal.The equation, expanded to also account for the effect of the NSB: •  where  is the mirror area, ( 0 ) is the quantum efficiency at the peak of the optical passband, ( 0 ) the optical efficiency of the rest of the system, ( 0 ) is the stellar differential photon flux, | | 2 ( 0 , ) is the squared visibility at the observed wavelength and baseline,   the effective cross-correlation electrical bandwidth,  the excess .Direct measurement of the angular stellar diameter of 9 reference stars, assuming a uniform disk profile, as a function of their reference diameter, also coming from a direct measurement over similar wavelengths.See Table 2 for their   ,   , assumed physical parameters and the source of the reference measurement.Dashed black line shows where the reference equals the measured diameter values.
noise factor of the PMTs,  the observation time,  the average background to starlight ratio and  the normalized spectral distribution of the optical passband (see eq. 5.6 in Hanbury Brown (1974)).More precisely, parameters , , ( 0 ), ,  and  must be understood as the geometric mean between the two telescopes.These parameters, first estimated by Acciari & al. (2020b), have been updated and are shown in Table 4.To test if the acquired data matches the sensitivity expected all data presented in this work was compared with the expected signal to noise inferred from eq. ( 6).As discussed in section 3, no quality cuts are applied to the data because a weighting is applied in the signal averaging step to account for the variable quality conditions of the observations.As shown in Fig. 11, once the signal weighting is taken into account when computing the expected signal to noise (by calculating a weighted observing time per correlation signal), the observed sensitivity matches the expectations.
We can also test the resulting relative uncertainty of the angular diameter measurements reported in this work, and compare them with the expected uncertainties.The resulting relative uncertainty of these measurements are shown as a function of stellar B magnitude in Fig. 12, before (grey circles) and after (coloured circles) correcting stellar B magnitude with the average atmospheric absorption.Note multiple effects may deviate our measurements from nominal performance: the different total exposure times acquired for each star, the different UV coverage acquired for different stars (MAGIC has only two telescopes at a fixed location) and the very different night-sky brightness during the observation of each star.However, the achieved relative uncertainty of the measured diameters are in good agreement with the expected trend.

Systematics evaluation
The systematic uncertainties associated with the analysis of data taken from VHE gamma-ray sources, i.e. the main scientific purpose of the MAGIC telescopes, have been studied in detail in previous works (Aleksić et al. 2012(Aleksić et al. , 2016b)).SII observations are not affected by most of the systematics generally associated to this technique: Uncertainties associated to the energy scale (number of photoelectrons detected), such as quantum efficiency, gain evolution effects due to the VCSELs, mispointing, untracked atmospheric transmission or light collection efficiency of the system only affect the S/N of the correlated signal, and not the measured Pearson's correlation.
Atmospheric turbulence is expected to add tiny differences in the actual path length of individual photons (Hanbury Brown 1974).This turbulence would ultimately become the limiting factor for   of an optical interferometer (in the case of perfect mirrors and photodetectors with ideal timing).The scale of these differences is expected to be below 30 ps, as calculated by Cavazzani et al. (2012), and therefore the effect is negligible for MAGIC-SII at the bandwidth scale of the current system.As expected, no significant broadening of MAGIC-SII correlation signals has been detected over high-zenith observations (up to 70 • in zenith distance).
Effects modifying the expected zero-baseline correlation (ZBC) of the system, or introducing a miscalibration in the visibility computation will enter the analysis as a systematic.The main sources of systematics that have been identified, summarised in Table 5, are the following: Stability of electronic bandwidth: As described in Hanbury Brown (1974), an untracked evolution of the electronic bandwidth of the system would lead to a modification of the ZBC, which needs to be stable for a reliable analysis.The single-photoelectron response of PMTs is considered stable, but as described in section 2.1, the signal transmission through VCSel's could perhaps add a time-dependent variation to the bandwidth.As described in section 2.4, the flexibility of the GPU correlator allows a simultaneous measurement of the cross-correlation between input signals and their auto-correlation, which translates into a simultaneous measurement of the electronic bandwidth of each channel.Auto-correlation measurements show the electronic bandwidth of the system is extremely stable over yearly timescales, with variations well below the 0.5% level (largest rare deviations on the 1% level, which can be easily identified).
Optical bandwidth: The optical bandpass, i.e. the resulting distribution in wavelength of the surviving photons after the narrow-pass optical filter in front of MAGIC PMTs, is the other parameter that would dominate a ZBC evolution.As described in section 2.2, the optical bandwidth of the system is dominated by the wide angle of incidence of photons reaching the filter.MAGIC AMC tracks the number and location of mirrors properly focused (which may vary between different observations).Dedicated toy MC simulations were employed to estimate the impact of missing mirrors over ( 0 ): 1) when missing mirrors are randomly distributed on the dish or 2) when removing 4% of the outer layer of mirrors the effect is still imperceptible; 3) when removing 4% of mirrors in the central part of the dish the effect is perceptible, but extremely small (< 1% level); 4) Only when removing a large fraction (more than 20 mirror facets) from the center of the dish the effect is larger than 1%, situation which never happened in the data presented here.Another effect contributing into the systematics, described by Hanbury Brown (1974), is the mirror deformation as a function of pointing elevation.Dish deformation biases the amount of mirrors actually illuminating the pixel, and may modify the true mirror area (and therefore photon incidence angle distribution on the narrow-band optical filter).In the case of the NSII, telescopes were equipped with optics collimating the beam, improving the performance of the optical filter but increasing the effect of this systematic.MAGIC, on the other hand, is equipped with AMC, which strongly reduces this effect, making it second order for full-mirror observations, as the location of unfocused mirrors is not dominated by the dish bending.
Gain evolution of DC readout: MAGIC slow control measures the DC of all camera pixels to track their illumination level and ensure their safety.As described in section 3.2, these DC measurements from MAGIC slow control are used to transform Pearson's correlation into visibility measurements.Therefore, any un-tracked gain variation within the DC ADC branch (PMT or DC ADC gains) will modify visibility measurements, and therefore will contribute into the systematics.The authors have identified several effects that could lead to gain variations in the DC branch: (i) temperature variations between summer and winter.Even if temperature within MAGIC cameras are controlled via cooling, the average temperature of the camera differs between seasons.Devoted lab measurements confirmed that the temperature compensation of the 12-bit ADC of the DC monitoring branch provides identical current measurements down to the part-per-million level.
(ii) PMT gains evolution over bright light source exposure.Lab test with pixel clusters identical to those used in the M2 camera indicate that when a PMT with a HV set in dark conditions is suddenly exposed to a bright light source (equivalent to a jump in current of ∼ 20 uA) the pixel gain undergoes a recovery period, asymptotically reducing the measured DC.From the controlled tests performed, this effect seems to be proportional to the charge applied, and is maximum when the PMT was previously kept in the dark, reaching significant gain variation peak amplitudes up to ∼ 8%, asymptotically reduced to 3%/2% after 5/10 minutes respectively.It should be noted that such gain variations can be tracked in our system and, under onsky observation conditions, this effect has never been observed with such an amplitude, and as it is only expected to affect significantly a small percentage of the observing time.We set a conservative overall limit over our measurements of ∼ 1% for very bright targets (used to compute  2 (0) ,  ).
(iii) PMT degradation.All PMTs undergo a slow process of gain degradation given the accumulated charge they are exposed to.We estimate, by assuming generous yearly charge exposures, a maximum yearly gain variation of 0.8%.Calibration measurements will be able to correct for this effect, but as best calibrators are mainly observable in winter, a conservative 0.8% systematic in the gain is assumed.
(iv) Observations with different HV will lead to different gains.A gain calibration is performed as described in section 3.2 to account for small HV variations.With a constant photon flux, () is expected to increase linearly with ().Very high current values are likely to escape this linear relation, and therefore would lead to different gains across our observations.To ensure this effect is not significant, calibration measurements over the HV and DC values we use for SII observations were performed, both to measure the relation between HV and DC for each pixel, and also ensure proportionality.No deviation from linearity was seen for any pixel over the safe current values employed.
Night-sky background DC subtraction: As described in section 3.2, the interferometry analysis needs to evaluate separately ( ) and ().To do so, we simultaneously gather DC measurements from our signal and background pixels.Conversion factors are needed to calculate the equivalent DC(NSB) at the signal pixel.There are several systematic sources associated to this background subtraction: a) If stars of a comparable brightness than the NSB enter background pixels FoV, their evaluation of the background will be overestimated; b) Slowly evolving effects may also affect this conversion factors, such as dirt and dust on the surface of individual narrow-band filters (not affecting the optical bandwidth, but would affect the relative DC measured between different pixels).As most of our observations (the full reflector focusing light to a single pixel) allow us to have one signal and two background pixels simultaneously, we can directly evaluate with the gathered data what is the scale of the systematic we expect from a wrong background evaluation.It should be noted that this is expected to be very small for bright stars, but it will increase for fainter stars, in which DC(NSB) will be a significant fraction of DC(Star).The impact on the evaluation of the DC(star) is expected to be 1.5/3% for bright/faint stars respectively (assuming faint stars are those fainter than ∼ 3.5 B mag).These are conservative values as during full-mirror observations we have 2 available background pixels: by selecting the one providing a lower DC(NSB), we suppress effect a), while by monitoring the ratio of DC(NSB) between background pixels allow us to update background-to-signal pixel conversion factors (and therefore, mitigate the effect b).
Residual electronic noise in the correlation: The MAGIC counting house, where the correlator is located, is equipped with multiple electronic devices producing high-frequency noise, that could eventually affect the interferometry setup.Any residual noise in the correlation would add systematics to the visibility measurements by deviating Pearson's correlation from the true value.We have done extensive tests to estimate the scale of this residual correlation, which can be seen when no photons are being recorded (HV off, and therefore signals exclusively come from electronic noise from digitizers and correlators).We evaluate the scale of this systematic by integrating over long observing periods, and evaluate the resulting deviations with respect to 0 for () far from the region where we expect the signal.When integrating very long observing times (> 50 h), a low-amplitude 1.67 MHz electronic noise is measured, stable in the () •  1 •  2 space (removing the normalization of Pearson's correlation).This systematic is strongly mitigated when adding a strong photon flux, and is therefore negligible for bright stars.This systematic, if not suppressed, would be significant when observing faint stars over dark conditions.By applying a 12 MHz cut in frequency (one order of magnitude narrower than our bandwidth) the residual electronic noise is strongly mitigated while correlation signals remain unaffected.
The stability of the zero-baseline correlation was evaluated by testing different procedures to measure  2 (0) ,  , each under different assumptions: i) Using eps CMa channel-pair wise data independently; ii) Using eps CMa data and assuming a common   ; iii) Using all available datasets from stars expected to be well described by a UD profile (removing fast rotators and spectroscopic binaries).In all cases,  2 (0) ,  are statistically consistent (largest variation of ∼ 2%).
It should be noted that the scale in time in which these systematics take effect may be very different, and therefore it will affect the reconstructed stellar diameter measurements in different ways.Long-term PMT degradation, even if corrected with calibration measurements updating  2 (0) ,  , will add a small systematic (< 0.8%) for summer sources (as most of the brightest reference stars are observable in winter).Effects expected to be transient over daily timescales (e.g.optical bandwidth deviations or some effects modifying PMT gains) are expected to cancel out when integrating datasets covering significantly longer time periods (e.g. when computing  2 (0) ,  , or when measuring faint sources with observations spanning over multiple weeks/months).For the measurement of  2 (0) ,  we use observations covering two different winters, each with multiple observing nights.In the case of eps CMa,  2 (0) ,  was measured simultaneously to , reducing the chances of systematic differences between  2 (0) and  2 () measurements, and therefore the only sources capable of systematically deviating the measured stellar diameter are gain drift and DC(NSB) subtraction effects (a maximum combined uncertainty of 1.8%).
The expected systematic uncertainty for each measured  were calculated by modifying squared visibility measurements in the most pessimistic way, and computing the limiting   , for each of these scenarios: • eps CMa: first/second half of visibility points were shifted up/down respectively by 1.8%.
• Others:  2 (0) ,  is shifted up/down by 1.8% while the rest of  2 measurements are shifted in the opposite direction by a brightnessdependent factor: 3.8% if   > 3.5 and 2% otherwise.
The reported systematic uncertainties in tables 2 and 3 refer to the difference between each   , and the best fit value, which is the largest deviation our systematics could cause in the most pessimistic scenario.

Future prospects
We have presented an evaluation of the performance of the current interferometry setup in MAGIC and the first measurements of diameters of several stars.As shown in Fig. 12, MAGIC can currently realistically target stars until ∼ 4 B mag, which means that we cannot compete yet with other long-baseline optical interferometers such as CHARA (Mourard et al. 2009).But these results prove the potential of planned improvements to boost sensitivity and angular resolution, described in more detail in Cortina et al. (2022).
Bandwidth improvements: As discussed in section 2.1, the current system bandwidth is limited by the digitizer used.We expect a factor ∼ √ 2 improvement in sensitivity by upgrading digitizers to the next generation .Improvements beyond this point would require faster photo-detectors and upgraded signal transmission.An alternative approach has also been tested during the last years: a faster readout based on a commercial 4 GHz sampling rate ADC coupled to an FPGA (Xilinx ZYNQ UltraScale+ RFSoC ZU28DR) which performs the correlation real-time.This readout currently operates in parallel to our nominal readout.Significant correlation signals have been detected although the system suffers from a strong correlated noise, with prominent components at frequencies beyond the bandpass of our nominal readout.
Photo-detection efficiency: a relatively straightforward way of improving the current system sensitivity would be to upgrade MAGIC photo-detectors to a newer generation of these detectors, as those used within the LST-1 camera.If such an array was implemented, the QE would increase from the current ∼ 32% to 40%.

Additional telescopes:
The first 23 m large-sized telescope prototype (dubbed LST-1) of the future Cherenkov Telescope Array Observatory (CTAO) was inaugurated in 2018 at a distance of ∼100 m from the MAGIC telescopes (Mazin & al. 2021).Three more telescopes with almost identical characteristics (LST2-4) are under construction and assembly, scheduled to become operational in 2025.
In addition to mirror area (370 m 2 ), photo-detectors QE (∼ 40%), and total optical efficiency (∼ 60%) are expected to be significantly better than the one from MAGIC.This increased photon-detection efficiency would dramatically improve sensitivity and the simultaneous UV coverage of the array.The expected relative errors as a function of stellar magnitude of MAGIC-LST1 and MAGIC-LST1-4 arrays is shown in Fig. 13 (Cortina et al. 2022).Minimal modifications in one cluster of LST-1 camera electronics have been designed and tested to enable optical transmission of a single pixel data to the MAGIC correlator.Adding additional telescopes to the array imply adding more input channels to the current GPU-based correlator.As described in section 2.4, the current system is already able to handle four input channels, which means no upgrades are necessary for the MAGIC-LST-1 array.However we typically operate with two pixels in all telescopes, to allow sub-reflector and chess-board observations.For this reason, the optimal implementation of adding the four LSTs would require a total of 12 input channels.A concept for a new GPUbased scalable correlator is currently being tested, described in detail by Cortina et al. (2022).If successful, this concept with be capable of handling all input signals coming from the full Northern CTAO.

CONCLUSIONS AND DISCUSSION
In this publication, we have demonstrated the feasibility of using IACTs as optical intensity interferometry arrays with only minor hardware additions.Given the narrow optical bandpass employed for these observations, they may be performed during full-Moon periods, times in which IACTs are generally not able to perform VHE gammaray astronomy.This breakthrough expands the scientific impact of IACT facilities and also increases their operational duty cycle.
We have found that the measurements made by the MAGIC-SII system are in good agreement with the measurements obtained by other observatories (see Table 2).Having demonstrated the capabilities of the system, we were thus able to reliably measure the diameters of 13 stars in the 400-440 nm band for the first time (see Table 3).Furthermore, we have shown that the sensitivity of the system is already capable of reaching relative errors at the few percent level over reasonable observing times (e.g.eps CMa or eta UMa).This level of precision, in principle, already allows for the study of the oblateness of fast-rotating stars.Measuring the oblateness of particular stars (such as blue supergiants, hypergiants, and Wolf-Rayet stars) allows to constrain their rotational speed, a crucial parameter that affects their structure and evolution.In order to allow for this or any other non-radially symmetric model fitting procedure, the MAGIC collaboration will release future MAGIC-SII data products using optical interferometry standards, such as OIFITS (Duvert et al. 2017).
In the near future, we will also be able to perform the first test observations by including the LST-1 in the array.The large increase in the expected S/N will realistically allow us to increase the precision of the diameter measurements, enabling us to achieve relative statistical errors below 1% over a large sample of stars.At that point, we expect to reach the scale of the systematics affecting our current analysis.The study of these systematics will also be easier by having a better S/N: shorter observation times will be required to make measurements of the correlation and, given the multiple simultaneous baselines and orientations covered, more stars can be employed to constrain them.
In conclusion, the results of the observations shown in this paper are consistent with expectations.Although they may not yet be competitive with other long-baseline interferometers (e.g., CHARA), this work demonstrates that the MAGIC-SII system will certainly be competitive if other telescopes belonging to the future Cherenkov Telescope Array Observatory northern array are included, as shown in Fig. 13.Adding additional telescopes to the array using the intensity interferometry technique is technically simpler than in the case of classical optical interferometry, which means it may become a realistic avenue for performing optical interferometry with intertelescope distances well beyond the km in the coming years.Although sensitivity may not extend beyond the 7-8 magnitude scale for such simple implementations (given the poor optical Point Spread Function of IACTs), the UV coverage achieved by large and dense arrays of IACTs would be unprecedented in the optical regime.

Figure 1 .
Figure 1.Diagram describing the different components of the MAGIC-SII setup.The main differences with respect to standard VHE observations are the incorporation of narrow-band optical filters positioned in front of MAGIC photodetectors and the continuous digitization and correlation of their recorded signals with the fastest bandwidth possible, performed by a set of digitizers and a GPU-based correlator.

Figure 2 .Figure 3 .
Figure 2. Electronic circuit design of the DC monitoring path.

Figure 4 .
Figure 4. Filter holder of the MAGIC telescopes below the white target, with room for six filters.At the time the picture was taken only two filters were installed (greenish circles), one of the holes was left open and three holes were closed with a plastic cap.Behind the filter holder and the diffusive white target one can see the hexagonal light concentrators (Winston Cones) right in front of the PMTs.

Figure 5 .
Figure5.Example of calibrated Pearson's correlation signal.Top) Weighted average correlation as a function of the zero-time corrected delay ( - 0 ) between the two input signals (A-C channels), from ∼ 30 min of observation with an average baseline of 52 m pointing to eps CMa (Adhara).The amplitude of the correlation is shown in the legend, using a Gaussian function with a 12 MHz cut high-pass filter applied.Best fit is drawn in blue while 1- uncertainty region is drawn in orange.Bottom) calibrated DC (in A) from each MAGIC telescope as a function of the zenith distance of the observation, both for the current associated to the star, as well as the NSB.Different DC values from the same telescope and zenith distance come from observations from different nights.Wind gusts produce relatively large fluctuations in M2 (less protected from wind than M1) while the low-amplitude sawtooth effect seen in both telescopes is due to the azimuth tracking.

Figure 6 .
Figure 6.Squared visibility vs baseline, showing the uniform-disc stellar diameter measurement of eps CMa using A-C correlation data (M1 pixel 251 correlated with M2 pixel 251).Envelopes surrounding each measurement highlight the true distribution of correlated baselines due to the large size of MAGIC reflecting dishes.Resulting best-fit values of the stellar diameter and normalization, i.e. the zero-baseline correlation  2 (0) − , are shown in the legend both from the  2 minimization as well as a bootstrapping procedure.Bottom figure shows the residuals with respect to the best fit uniform disc model.

Figure 7 .
Figure 7. Correlation signal delay (left) and width as its Gaussian  (right) for all V 2 measurements presented in this work.The dashed vertical line is the expected width of the correlation discussed in section 3.2.

Figure 8 .
Figure 8. Channel-pair-wise stellar diameter measurements relative to the combined reconstructed value.Only measurements of stars having at least two valid measurements with different channel pairs are shown.
Figure 9. Direct measurement of the angular stellar diameter of 9 reference stars, assuming a uniform disk profile, as a function of their reference diameter, also coming from a direct measurement over similar wavelengths.See Table2for their   ,   , assumed physical parameters and the source of the reference measurement.Dashed black line shows where the reference equals the measured diameter values.

Figure 10 .
Figure 10.Direct measurement of angular stellar diameter of 13 stars newly measured by MAGIC, assuming a uniform disk profile, as a function of their expected angular diameter from Bourges et al. (2017).See Table3for their   ,   and assumed physical parameters.Note the large uncertainty associated to phi Sgr is due to the very small observation time acquired (15 min).

Figure 11 .
Figure 11.Comparison between the measured and expected signal to noise of each correlation signal.Grey points use the total observing time of each observation to compute the expected signal to noise, while green points take into account the weighting applied (described in section 3) to minimize the impact of low-quality data (observations with significantly lower photon flux).

Figure 12 .
Figure 12.Relative uncertainty of measured stellar diameters as a function of their B magnitude.The diameter of each point is proportional to the measured stellar diameter, and the colour scale shows the total observing time used (no quality cuts applied).The nominal expected error as calculated by Cortina et al. (2022) is shown as a dashed line, and as calculated by Fiori et al. (2022) as a dashed-doted line.Grey points use the average B magnitude, while coloured ones are corrected for the atmospheric extinction affecting the observations.

Figure 13 .
Figure 13.Relative uncertainty of measured stellar diameters as a function of their B magnitude for the current MAGIC interferometer, an extension to MAGIC and LST-1, and a further extension to MAGIC and the four LSTs in La Palma.Dashed lines indicate sensitivity during dark time, while solid points estimate the sensitivity under conservative very-high NSB illumination levels (both for a 2.5 h observing time).The simulated stellar position and diameter is the one from gam CrV.

Table 1 .
Response of individual MAGIC pixels to calibration pulses as measured by the MAGIC-SII GPU-based correlator.Telescope, pixel ID number, input and measured average widths (Gaussian ) for each pixel currently equipped for interferometry observations.These are computed by averaging calibration pulses from raw digitized data (500 Msa/s).

Table 2 .
Table of reference stars sorted by B magnitude.Spectral types and B magnitudes are from SIMBAD.The reference angular diameters are from the NSII (Hanbury

Table 4 .
Estimated performance parameters of the MAGIC-SII setup.

Table 5 .
Evaluated systematic uncertainties over squared visibility measurements identified to effect the MAGIC-SII system.