Radiative Properties of Plasmoids and Plasmoid Mergers in Magnetic Reconnection

Blazars are relativistic collimated plasma outflows, also known as jets, that point very close to our line of sight. These jets are launched from accreting supermassive black holes with very high magnetic energy (Blandford & Znajek, 1977; Tchekhovskoy et al., 2010). Around sub-parsec distance to a few parsecs, there is an unresolved region that emits highly variable nonthermal-dominated multi-wavelength emission, often referred to as the blazar zone. The variability time scale can reach as short as a few minutes (Albert et al., 2007; Aharonian et al., 2007; Ackermann et al., 2016), indicating extremely fast particle acceleration. Multi-wavelength campaigns often find that the spectral index of some observational bands hardens when those bands exhibit flares (Giommi et al., 1990; Abdo et al., 2010; Krauß et al., 2016), referred to as the “harder-when-brighter” trend. This trend further supports that these blazar flares are indeed driven by particle acceleration in situ in the blazar zone. The radio to optical blazar emission, in some cases up to X-ray bands, is dominated by the synchrotron emission from ultra-relativistic electrons, evident by the observed high PD (Scarpa & Falomo, 1997; Smith, 2017; Lister et al., 2018). Observations find that the polarization signatures can also evolve in time (Chandra et al., 2015; Covino et al., 2015; Itoh et al., 2016; Jorstad et al., 2022). Very interestingly, optical polarization monitoring programs find that PA swings, defined as angle rotating continuously in one direction for $\gtrsim 90^{\circ}$, are typically accompanied by multi-wavelength flux variability (Marscher et al., 2010; Larionov et al., 2013; Blinov et al., 2016, 2018), hinting the co-evolution of magnetic fields during particle acceleration (Zhang et al., 2015; Zhang et al., 2018). Such PA swings have been reported in soft X-rays by the recently launched IXPE as well (Di Gesu et al., 2023). It is therefore believed that the blazar zone is probably magnetized, where the jet dissipates a significant amount of its magnetic energy to accelerate particles. Although the radio to optical emission is dominated by nonthermal electrons, the blazar jet may contain a proton population, which is evident by the detection of a very high energy neutrino event coinciding with the blazar TXS 0506+056 flare (IceCube Collaboration et al., 2018). Understanding the physical conditions of the blazar zone and particle acceleration mechanisms therein is of critical importance to the time-domain and multi-messenger astronomy.

Magnetic reconnection is the primary plasma physical mechanism to accelerate nonthermal particles in a magnetized environment. This plasma physical process happens when two oppositely oriented magnetic field lines come close to each other, break and rejoin to rearrange the magnetic topology, and dissipate a large amount of the magnetic energy between the field lines. Such magnetic field topology can form during magnetic instabilities such as kink instabilities or in a striped jet configuration (Begelman, 1998; Giannios & Spruit, 2006; Giannios & Uzdensky, 2019; Alves et al., 2018; Davelaar et al., 2020). Previous particle-in-cell (PIC) simulations show that reconnection can efficiently accelerate both electrons and protons (Sironi & Spitkovsky, 2014; Sironi et al., 2015; Guo et al., 2014; Guo et al., 2016, 2021; Werner et al., 2016, 2018; Li et al., 2018; Cerutti et al., 2012; Zhang et al., 2021b; Chernoglazov et al., 2023), making it an ideal particle acceleration mechanism for hadronic blazar models. Particles are accelerated via Fermi-like mechanisms (Guo et al., 2023; French et al., 2023), or they may be accelerated directly in the electric field generated by reconnection (Sironi, 2022). The accelerated particle spectral indices generally depend on the magnetization factor, which is the ratio of magnetic energy over enthalpy, $\sigma$, and the guide field strength, which is the magnetic field component that is perpendicular to the anti-parallel field lines, $B_{g}$ (Guo et al., 2015; Sironi et al., 2015). Combined PIC and radiation transfer simulations find that reconnection typically leads to the harder-when-brighter trend during flares and polarization can strongly vary as well (Zhang et al., 2020; Zhang et al., 2021a). When the 2D reconnection starts, the current sheet, the thin layer between the anti-parallel magnetic field lines, can generate a series of quasi-circular structures with concentrated magnetic fields and nonthermal particles, called the plasmoids (Sironi & Spitkovsky, 2014; Guo et al., 2014; Petropoulou et al., 2018; Hakobyan et al., 2021). They contribute the majority of nonthermal emission (Zhang et al., 2018; Christie et al., 2019). Plasmoids move in the outflow direction of the reconnection region. Some plasmoids can reach relativistic speeds; if the outflow direction is along the line of sight, they can lead to extra relativistic Doppler boosting (the jet-in-jet model) and produce ultra-fast flares (Giannios et al., 2009). Plasmoids can also collide into each other and merge, during which they can dissipate a large amount of their magnetic energy and accelerate additional nonthermal particles, resulting in flares and PA swings (Zhang et al., 2020, 2022). Although 3D simulations suggest that reconnection can generate strong turbulence that results in more complicated structures, 2D slices along the guide field direction still generally show plasmoid structures (Guo et al., 2021; Huang & Bhattacharjee, 2016). Given the very high computational cost, it is impossible to simulate a realistic blazar zone using combined PIC and radiation transfer to study the radiation signatures of magnetic reconnection. It is therefore important to understand the radiative properties of plasmoids in more detail to interpret the radiation signatures on realistic physical scales.

Sironi et al. (2016) is a pioneering study on the physical properties of plasmoids. Using the local saddle points of the magnetic vector potential, the paper segments plasmoids and finds that plasmoids are generally self-similar structures. This paves the way for modeling radiation signatures of magnetic reconnection on the scale of realistic blazar zones. Christie et al. (2019) models the multi-wavelength blazar radiation from reconnection based on the above findings. However, the above segmentation method has a major drawback: by definition, adjacent plasmoids or groups of plasmoids share a common saddle point. Thus all plasmoids appear to be merging, making it impossible to distinguish isolated and merging plasmoids and study their respective radiative properties. Moreover, the radiation signatures in Christie et al. (2019) are not computed directly from PIC simulations, but based on simplified modeling assumptions inspired by Sironi et al. (2016). For instance, one of the key assumptions is to consider plasmoids as unresolved concentrations of magnetic fields and particles, preventing the study of polarization signatures. To further study the radiative properties under first principles and enable direct comparison with 3D simulations and in future large-scale hybrid simulations, we need to better segment plasmoids and their mergers to explore their radiative properties statistically.

Here we introduce a new segmentation method using image processing techniques to identify isolated plasmoids and plasmoid mergers in combined 2D PIC and radiation transfer simulations of magnetic reconnection. The boundaries of plasmoids are determined by both the magnetic vector potential and particle density. Our approach can distinguish radiation contributions from plasmoids, mergers, and the primary current sheet. By statistically studying hundreds of plasmoids and tens of mergers from nine simulations with different physical parameters, we are able to examine their radiative properties, including polarization signatures. Section 2 presents our simulation setup and methodology, Section 3 shows the overall radiation signatures based on our segmentation method, Sections 4 and 5 describe radiative properties of isolated and merging plasmoids, respectively. We summarize and discuss our results in Section 6.

The goal of this paper is to study the radiative properties of isolated and merging plasmoids statistically. Thus we include radiative cooling effects for synchrotron and Compton scattering in our simulations (Zhang et al., 2022). Sironi et al. (2016) shows that the open boundary simulation is advantageous in generating a large number of plasmoids, but the number of plasmoid mergers, especially mergers between large plasmoids, is very limited. Here we choose periodic boundary conditions for our combined PIC and radiation transfer simulations, which guarantee that most of plasmoids will merge so as to boost the plasmoid merger statistics. We use image processing techniques to combine both the magnetic vector potential and plasma density to identify the boundary of plasmoids and their mergers. We find that our approach can clearly separate isolated plasmoids, plasmoids mergers, and the primary current sheet. Additionally, both isolated and merging plasmoid boundaries can nicely envelop their synchrotron radiation, ideal for studying the radiative properties. Our method can also track the evolution of both isolated and merging plasmoids. We are able to segment nearly one thousand plasmoids and a few hundred mergers from nine simulations with various physical conditions. We find that the PA swing driven by mergers is due to the merging plasmoids rotating around each other, resulting in the newly accelerated particles at the merging site streaming along the curved magnetic field lines surrounding the merging plasmoids as found in previous works (Zhang et al., 2018, 2021a). This section presents our simulation setups and segmentation methods. We also describe how we define the radiative properties of isolated and merging plasmoids in the statistical studies.

This section shows the overall radiation signatures from isolated and merging plasmoids as well as the primary current sheet. We find that the emission contribution from the primary current sheet is trivial compared to the isolated and merging plasmoids. Stokes $I$, $Q$, $U$ during flares are dominated by merging plasmoids, although this is probably biased due to the periodic boundary conditions, which can lead to a lot of mergers. All continuous PA rotations $\gtrsim 90^{\circ}$ are driven by mergers, and for swings $\gtrsim 180^{\circ}$, the Stokes $Q$ and $U$ follow roughly sine functions with $\sim 90^{\circ}$ difference in phase. Isolated plasmoids can have varying PD and PA due to the non-uniform distribution of nonthermal particles inside, but they cannot create $\gtrsim 90^{\circ}$ rotations.

Figure 5 shows the overall temporal evolution of the optical Stokes $I$, $Q$, $U$ as well as PD and PA. The general patterns are similar to Zhang et al. (2020): since the reconnection starts, the optical light curve (Stokes $I$) shows multiple flares till the reconnection gradually saturates at $t\gtrsim 2.5\tau_{lc}$. The first four relatively short flares before $t\sim 1.5\tau_{lc}$ are mostly due to isolated plasmoids and mergers in the reconnection layer, while the last giant flare between $t\sim 1.75\tau_{lc}$ and $t\sim 2.75\tau_{lc}$ results from a large binary merger. The first PA swing near $t=1.2\tau_{lc}$ is due to mergers between multiple plasmoids, while the sudden PA drop around $t=2.25\tau_{lc}$ is due to the plasmoid merger shown in the previous section. Previous works have shown significant changes in PD and PA during flares (Zhang et al., 2018; Zhang et al., 2020; Zhang et al., 2021a), but here we also present how the Stokes $Q$ and $U$ vary to further explore the physical origin of polarization variations. As we can see, major rises and falls of the PD are highly correlated to flares of Stokes $Q$. The Stokes $U$, on the other hand, do not show a clear correlation with the PD or PA. The reason lies in the radiative properties of isolated plasmoids and mergers. Isolated plasmoids, due to their geometric shape, have a net polarization that is perpendicular to the reconnecting magnetic field lines (described in Section 4). On the other hand, the secondary reconnection at the contact region of the merging plasmoids, generally leads to polarization along the reconnecting field lines. These two polarization directions give the positive and negative Stokes $Q$, while the Stokes $U$ results from particles accelerated at the secondary reconnection region that stream along the post-merger plasmoid, which have been partially cooled via radiation.

Our segmentation method allows us to distinguish different emission contributions. Figure 5 shows that in Stokes $I$, $Q$, and $U$, plasmoid mergers contribute the most to emission, isolated plasmoids contribute a small portion, while the current sheet has a minimal contribution. We note that the dominance of plasmoid mergers in all Stokes parameters in our simulation is due to the periodic boundary conditions, which result in many more mergers compared to open boundaries as in Sironi et al. (2016). Nonetheless, at an earlier stage ($0.5\tau_{lc}<t<1.5\tau_{lc}$), when both isolated and merging plasmoids exist in the simulation domain, merging plasmoids still contribute much more than isolated ones in radiation. This is because plasmoid mergers not only consume a significant amount of the magnetic energy in participating plasmoids, they also harden the particle spectra (described in Section 5), considerably boosting the optical flux. Radiation from the current sheet appears to be significant at the beginning of the reconnection ($t\lesssim 0.75\tau_{lc}$); however, this is because we do not include plasmoids smaller than $16\times 16$ PIC cells, whose radiation is counted as from the current sheet. If we consider any plasmoids that are resolved by at least 2 PIC cells, then even the emission in the beginning of the simulation is dominated by isolated plasmoids (about 70% to 80%) and mergers (the rest), while the contribution from the current sheet is trivial. The Stokes $Q$ and $U$ are dominated by the mergers as well, further supporting that the polarization variations result from plasmoid mergers. We note that although mergers almost continuously happen during our simulations, there are a few snapshots where very few mergers are ongoing. In those snapshots, radiation from isolated plasmoids becomes significant, which leads to sudden changes in their contribution to the Stokes $I$, $Q$, $U$ as shown in Figure 5. Such sudden changes are simply due to the periodic boundaries that make mergers so frequent: if the simulation has open boundaries, where plasmoid mergers are much less frequent, the transition will appear smooth.

Figure 6 presents a different run, in which a smooth large-amplitude PA swing happens. The general patterns are similar to Run 1, but this run does not have a major plasmoid merger near the end of the simulation. Instead, all plasmoid mergers are complete before $t\sim 1.75\tau_{lc}$, leaving two giant plasmoids that are far apart from each other (however, if the simulation continues for a longer time, they will merge eventually). As a result, the emission is dominated by isolated plasmoids later in the simulation. The distributions of nonthermal particles and magnetic fields can evolve in both plasmoids. They produce fluctuations in Stokes $Q$ and $U$, which lead to changes in the PA. But such fluctuations are not significant to make any PA swings $\gtrsim 90^{\circ}$. Since the two isolated plasmoids are the only dominating emission patches, the PD appears very high when Stokes $Q$ and $U$ rise. However, the reconnection has mostly stopped at this stage and the flux is very low, such high PDs are unlikely observable. The main difference of this run from Run 1 is the smooth large-amplitude PA swing around $t\sim 1.5\tau_{lc}$, which is due to a major plasmoid merger. We can see that the Stokes $Q$ follows roughly a cosine curve while the Stokes $U$ follows roughly a sine curve. The $Q$ and $U$ thus complete a quasi-circle around zero in the Stokes $Q$-$U$ plane (not shown here), a clear signature for PA swing found in observations. We will quantify the radiative properties of plasmoid mergers in Section 5. As we will see, the PA rotation amplitudes for plasmoid mergers can be stochastic, but they follow a general trend that depends on the size ratio of the merging plasmoids.

Here we present the physical and radiative properties of isolated plasmoids. We find that the majority of physical properties of isolated plasmoids are similar to those found in Sironi et al. (2016). The difference is that our simulations include radiative cooling, thus the kinetic energy of particles in plasmoids can decrease, resulting in a lower equipartition factor towards the end of simulations (refer to Equation 1). Similarly, the spectral peak of the nonthermal electrons decreases due to cooling (refer to 7). We find that all isolated plasmoids, regardless of their sizes, share approximately the same emissivity. They also exhibit a nonzero PD due to their quasi-circular geometric shape. We note that due to the periodic boundaries in our simulations, the number of isolated plasmoids is much less than that in Sironi et al. (2016). Thus we use Runs 1, 4, 5, 6, and 7 that share the same magnetization factor $\sigma$ and cooling factor $C_{10^{4}}$ to increase the statistics. Unfortunately, our image processing techniques are still too computationally expensive to explore dependence on additional parameters. Nonetheless, key physical and radiative properties of plasmoids do not show clear dependence on their sizes, indicating that isolated plasmoids are indeed self-similar structures.

Figure 7 shows the 2D histogram of four physical properties (in the comoving frame of the simulation domain) of isolated plasmoids against the plasmoid widths $w$ in the $z$ direction (i.e., perpendicular to the primary current sheet). As shown in the upper left panel, the plasmoids generally appear quasi-circular, where the width $w_{\parallel}$ along the primary current sheet (i.e., the $x$ axis) slightly larger than the width $w$ perpendicular to the primary current sheet, in agreement with Sironi et al. (2016). The stretching along the primary current sheet is quite obvious in Figure 1 as well. The ratio $w_{\parallel}/w\sim 1.6$ (red dashed line in the upper left panel) is more pronounced for larger plasmoids that are well resolved. Small plasmoids of $w\lesssim 0.02L$ ($w=0.02L$ is marked by the white vertical dotted line in Figure 7) have a radius that contains only $\sim 20$ PIC cells. In this case, the systematic errors in the plasmoid contours can become non-trivial, thus there are more plasmoids $w_{\parallel}/w$ deviating from the 1.6 ratio.

The lower left panel of Figure 7 shows the ratio of the average electron density in the plasmoids $n_{e}$ over that in the plasma inflow $n_{0}$. Our result is slightly higher than that in Sironi et al. (2016): our total magnetization factor is $\sigma\sim 22$, but we find $n_{e}/n_{0}\sim 5$ (red dashed line in the lower left panel), which is the value for the $\sigma=50$ simulation in Sironi et al. (2016). This is likely due to our different setups: we consider a proton-electron plasma with guide fields and radiative cooling, thus it is not an apples-to-apples comparison with Sironi et al. (2016). Our value seems to match well with the $n_{e}/n_{0}\propto\sqrt{\sigma}$ scaling in Lyubarsky (2005); however, we are unable to verify the trend with more runs with different $\sigma$ due to the very computationally expensive image processing method. We note that our results are still consistent with the $n_{e}/n_{0}\gtrsim\Gamma$ argument in Sironi et al. (2016), but we remind the readers that due to the periodic boundary conditions, our plasmoids are at most mildly relativistic, $\Gamma\lesssim 2$.

Other physical quantities (not shown) are in good agreement with Sironi et al. (2016), just like $w_{\parallel}/w$ and $n_{e}/n_{0}$. However, our simulations include radiative cooling, hence the average electron kinetic energy in plasmoids can vary in time. This is best represented by the equipartition factor $\chi_{eq}$ shown in the upper right panel. We define $\chi_{eq}$ in a similar way as in Sironi et al. (2016),

but we include both the kinetic energy density of protons and electrons in $\varepsilon_{\rm kin}$ since we use proton-electron plasma, while $\varepsilon_{\rm B}$ denotes the magnetic energy density. As stated in Sironi et al. (2016), the strong kinetic dominance for small plasmoids can be attributed to that field lines generally wrap around large plasmoids, so that the small plasmoids preferentially lie in regions where the field lines are less dense, resulting in weaker magnetic fields. However, large plasmoids in our simulations reach a state that the magnetic energy is slightly dominating ($\chi_{eq}\sim 0.4$), different from $\chi_{eq}\sim 0.6$ in Sironi et al. (2016). This is mainly due to radiative cooling. Previous works have found that in proton-electron plasma, the energy in nonthermal electrons and protons are comparable, i.e., $\varepsilon_{kin,e}\sim\varepsilon_{kin,p}$, and plasmoids will reach equipartition ($\varepsilon_{kin,e}+\varepsilon_{kin,p}\sim\varepsilon_{B}$) if there is no cooling. Since large plasmoids take a long time to grow, nonthermal electrons therein can radiatively cool and lose kinetic energy. Protons, due to their much longer cooling time, experience trivial cooling effects. As a simple estimate, if the electrons lose all their kinetic energy due to cooling, then the equipartition factor will be approximately $\chi_{eq}\sim\varepsilon_{kin,p}/(\varepsilon_{kin,p}+\varepsilon_{B})\sim 1/3$. In practise, electrons will only lose part of their kinetic energy, thus the $\chi_{eq}$ for large plasmoids is higher, at roughly 0.4 (red dashed line) in Figure 7 upper right panel.

The electron cooling is more clearly illustrated in the lower right panel of Figure 7, where we plot the spectral peak of the electron spectrum (peak refers to the peak in a log-log plot of $\gamma_{e}^{2}N_{e}(\gamma_{e})$ against $\gamma_{e}$). We find that small plasmoids can reach $\gamma_{peak}\sim 10^{4}$, which is very close to the electron magnetization factor $\sigma_{e}=4\times 10^{4}$. However, this value drops nearly linearly for larger plasmoids (red dashed line in the lower right panel). Since $\gamma_{peak}$ cools linearly with time, this trend infers that the plasmoid size grows linearly with time, i.e., a power law of index -1 in Figure 7, consistent with Sironi et al. (2016). Note that the plasmoid growth in our simulations consists of both the growth as isolated plasmoids and merging with other plasmoids, indicating that both are linear in time.

The luminosity of plasmoids generally scales with the square of the plasmoid width, as shown in the red dashed line in Figure 8 upper panel. This suggests that the emissivity in plasmoids is roughly constant regardless of the plasmoid size. We note that the luminosity of small plasmoids (smaller than the white dotted line) tends to deviate to slightly lower luminosity than the red dashed line, even though they have the highest $\gamma_{peak}$. This is because the small plasmoids have relatively weaker magnetic fields since they preferentially stay in regions with less dense magnetic field lines as mentioned above. As a result, their luminosity is slightly lower.

Due to the quasi-circular shape of the plasmoids, they yield a nontrivial PD. In our simulations, the viewing angle is along $+y$, the expected PD can be estimated via the following integration in an ellipse. For simplicity, we assume that the magnetic field strengths and nonthermal electron distributions are constant in the quasi-circular plasmoid, then the integration of Stokes $I$, $Q$, $U$ are given by

where $s=\epsilon$ for Stokes $I$, $s=\epsilon\Pi\cos{2\alpha}$ for Stokes $Q$, and $s=\epsilon\Pi\sin{2\alpha}$ for Stokes $U$, and $\Pi=70\%$ is assumed to be a constant for the synchrotron PD at every point within the ellipse as the electron power-law distribution has an index of $p=-2$, while $\alpha$ is the angle between the line perpendicular to the tangent of the ellipse at $\theta$ and $x$-axis. Since $s$ does not depend on $r$, we have

where $r=ab/\sqrt{(b\cos{\theta})^{2}+(a\sin{\theta})^{2}}$ is the radius of the ellipse from the center at different $\theta$ in which $a=w_{\parallel}/2$ and $b=w/2$. Therefore, we have

With some math one can find $\alpha$ as

Then Stokes $Q$ and $U$ are given by

By observing the expression for Stokes $U$, the integral is proportional to $\tan\theta$, thus the integration is zero. The Stokes $Q$ integration can be done analytically by noticing that $d\tan\theta=\sec^{2}\theta d\theta$, and we use $z=\tan\theta$, then we have

Thus we have

For $w_{\parallel}=1.6w$, we find that PD is about 16%. If we view from a different angle, the Stokes $Q$ and $U$ in the above equation need to be corrected by a projection factor. The net PD is zero only if we are viewing in the $xy$ plane along an angle $\sim\arctan{1/1.6}$ from $x$ axis, where the projected ellipse becomes a perfect circle. In practice, the nonthermal particles and magnetic field distributions inside isolated plasmoids can change in time, causing fluctuations in PD and PA as shown in the previous section. The above PD estimate should be considered as the average value. The PD deviates from the theoretical value for small plasmoids as shown in Figure 8 lower right panel. This is because our radiation transfer code has a lower resolution than PIC. For plasmoids smaller than $w=0.02L$ (the white vertical line), the plasmoid $w$ and $w_{\parallel}$ are only resolved by $\sim 3$ radiation cells, less than 9 cells for the plasmoid area. This should yield a PD of $\sim 1/\sqrt{9}\times 70\%=23.3\%$ even if the magnetic fields are completely random in the plasmoids. The smallest plasmoids contain only one radiation cell, resulting in the maximal PD $\sim 70\%$.

This section presents the radiative properties of binary plasmoid mergers. We find that the flare amplitude is roughly anti-proportional to the size ratio of the merging plasmoids. Additionally, the optical photon spectral index becomes harder by $\sim 0.5$ during the merger if the width ratio of the merging plasmoids is $R_{w}\lesssim 5$. Similarly, PA variations only occur with the same width ratio, whose amplitudes can lie between 0 and $180^{\circ}$ with an average of $90^{\circ}$. Finally, the merging time scale is proportional to the size of the secondary current sheet, with a rate of $\sim 0.1c$.

As shown in the previous section, large plasmoids have more energy than small plasmoids. As a result, mergers between relatively large plasmoids generally dissipate more energy and radiate more strongly than mergers between relatively small plasmoids. Here we choose to study the relative changes in the flux, optical photon spectral index, and PA to avoid differences due to the size of plasmoids participating in the merger. Additionally, we study the dependence of these changes on the width ratio of the merging plasmoids instead of their widths. Figure 9 (left panels) shows the results¹¹1Several mergers with width ratio $R_{w}>13$ are not shown as they lead to no changes in flux or polarization.. With the above choices, we find that indeed as long as the width ratio is similar, mergers leading to large post-merger plasmoids, which happen between relatively large plasmoids, yield similar variations in radiation signatures compared to those making small post-merger plasmoids (except for $\Delta\alpha$, where the small plasmoid mergers show smaller spectral changes than large ones). Additionally, Runs 2 and 3 with different $\sigma$ and Runs 8 and 9 with different cooling factors do not show clear difference from the other runs, since we only consider the relative changes in the flux and polarization during plasmoid mergers. Therefore, we only consider the width ratio parameter here and include all runs to boost statistics in the number of mergers. To better illustrate the statistical results, we bin the mergers based on the width ratio, and the right panels of Figure 9 show a skewed Gaussian distribution for the statistics: the 0.84, 0.5, and 0.16 quantiles of the skewed Gaussian distribution are plotted as the upper y error bar, median, and lower y error bar, respectively.

We find that only mergers with a width ratio $R_{w}\lesssim 3$ can produce a flare with an amplitude greater than 1 (i.e., double the flux). Mergers with a large width ratio $R_{w}\gtrsim 8$ bring about trivial flare regardless of the size of participating plasmoids. Interestingly, the flare amplitude seems to follow an anti-proportional relation to the width ratio. This can be understood in the following way. As shown in the previous section, isolated plasmoids have approximately the same emissivity. For a binary merger, say the larger participating 2D plasmoid has a surface area of $A$, then that of the smaller plasmoid is $\frac{1}{R_{w}^{2}}A$, where $R_{w}$ is the width ratio of the larger plasmoid over the smaller one. Therefore, the luminosity of the two isolated plasmoids right before the merger starts is given by

When the two plasmoids start to merge, the length of the contact region obviously cannot be larger than the width of the smaller plasmoid, thus the size of the secondary reconnection layer is at most proportional to the size of the smaller plasmoid. Since the magnetic energy density is approximately the same for plasmoids of any sizes, if a constant portion of the magnetic energy in the secondary reconnection layer is dissipated and transformed to radiation, then the extra luminosity induced by the merger should be $L_{1}=C\frac{1}{R_{w}^{2}}A$, where $C$ describes a portion of the magnetic energy of the smaller plasmoid that is dissipated. As a result, the flare amplitude is given by

This fitting curve is shown in red in the upper right panel of Figure 9. It is clear that this model describes the flare amplitude very well. However, it deviates slightly from the statistics for mergers of width ratio $R_{w}\sim 1$. This is because for $R_{w}\sim 1$, the secondary reconnection region is no longer proportional to the size of the smaller plasmoid, but it has to be smaller than that. As a result, $L_{1}<C\frac{1}{R_{w}^{2}}A$, so that the actual flare level is lower than the theoretical value of the red curve in Figure 9.

The optical photon spectral index hardens by $\sim 0.5$ for $R_{w}\lesssim 5$. Along with the flare amplitude, plasmoid mergers with $R_{w}\lesssim 3$ can well explain the often observed harder-when-brighter trend. This hardening feature quickly disappears for $R_{w}>5$. Obviously, if the flare level is too low, generally we do not expect to see spectral hardening. But to understand this behavior in more detail, we compare the electron spectrum at the beginning of the merger and at the flare peak of the merger. The main difference is that at the flare peak, the electron spectrum nearly extends to the electron magnetization factor $\sigma_{e}$ with an index of $p\sim 2$, as shown in the upper panel of Figure 4 (blue is near the beginning while the red curve is near the peak). This suggests that regardless of the width ratio of the participating plasmoids, mergers result in a hard spectrum of newly accelerated particles. However, if the width ratio is too large, the total number of newly accelerated particles is smaller than the number of particles already existing in the plasmoids participating in the merger, then the total particle spectral shape is trivially affected by the newly accelerated particles during the merger. This explains why $\Delta\alpha$ quickly drops to zero for width ratio $R_{w}>5$. For $R_{w}\lesssim 5$, $\Delta\alpha$ has a relatively large range because of the spectrum in the plasmoids. As shown in the previous section, small plasmoids, which are newly generated, have an electron spectrum extend to $\gamma_{e}\sim 10^{4}$, similar to the particle spectrum accelerated during mergers. However, large plasmoids, due to cooling, have lower spectral cutoffs. Since the optical synchrotron emission is dominated by electrons at $\gamma_{e}\sim 10^{4}$ in our setup, large plasmoids have the optical emission deep in the radiatively cooled spectrum. Consequently, when they merge, the newly accelerated particles, which are not yet cooled, push the optical photon spectrum much harder. Small plasmoids, on the other hand, have hard optical photon spectra before the merger anyway, thus the merger does not harden the spectral index significantly. This explains why mergers between small plasmoids tend to occupy a region with smaller $\Delta\alpha$ even for $R_{w}\lesssim 5$ in Figure 9 middle left panel.

PA variations are only considerable for $R_{w}\lesssim 5$ for similar reasons: PA changes are induced by newly accelerated particles; if they are not enough to affect the nontrivial PD resulting from the quasi-circular shape of the plasmoids, PA will not change significantly. Most of the PA variations are $\sim 90^{\circ}$. The reason is that the plasmoid mergers typically start from head-on collision between two plasmoids, so that the newly accelerated particles are at the secondary reconnection layer that is perpendicular to the primary current sheet. The magnetic field lines in the secondary current sheet are then parallel to the net polarization direction of the isolated plasmoids, thus the newly accelerated particles preferentially emit synchrotron polarized in the $\sim 90^{\circ}$ from the polarization direction of the pre-merger plasmoids. To complete a full $\sim 180^{\circ}$ PA rotation, the merger needs to meet three criteria: the merger cannot complete before the two plasmoids rotate $180^{\circ}$ from each other; the newly accelerated particles at the secondary current sheet cannot be mostly cooled before the rotation completes; the secondary current sheet cannot generate too disordered magnetic fields. For instance, the merger shown in Section 2 fails to make the $\sim 180^{\circ}$ PA rotation because of the very disordered magnetic field in the secondary reconnection region. Therefore, flares driven by reconnection are not necessarily accompanied by large PA swings.

Figure 10 shows the growth rate due to the merger. This is simply the width of the post-merger plasmoid minus the width of the larger plasmoid that participates in the merger, plotted against the duration of the merger defined in Section 2. Due to the limited number of mergers between large plasmoids, we do not have very good statistics for the growth rate. But Figure 10 can be fit with a linear growth rate that is $\sim 0.1c$, consistent with the growth rate of isolated plasmoids found in Sironi et al. (2016).

To summarize, we statistically study the radiative properties of isolated and merging plasmoids in magnetic reconnection. Our results are based on combined 2D PIC and polarized radiation transfer simulations. Our image processing method that combines the magnetic vector potential and particle density can clearly identify the boundaries of isolated and merging plasmoids. Our key findings are as follows.

1. 1.

   Radiation from magnetic reconnection is dominated by isolated and merging plasmoids. The contribution from the primary current sheet is mostly negligible.

2. 2.

   The luminosity of isolated plasmoids is proportional to its size.

3. 3.

   Isolated plasmoids have comparable magnetic energy density and particle kinetic energy, but the former is slightly higher.

4. 4.

   Isolated plasmoids generally have non-trivial polarization due to their quasi-circular morphology.

5. 5.

   Flares driven by plasmoid mergers have amplitudes that are anti-proportional to the size ratio of the plasmoids participating in the mergers.

6. 6.

   Mergers between plasmoids with a relatively small width ratio ($R_{w}\lesssim 5$) can lead to spectral hardening near the synchrotron spectral peak.

7. 7.

   Mergers between plasmoids with relatively small width ratio ($R_{w}\lesssim 5$) can produce PA variations ranging from 0 to $\sim 180^{\circ}$, with an average of $\sim 90^{\circ}$.

8. 8.

   PA swings are due to nonthermal particles that stream along the magnetic field lines surrounding the merging plasmoids, which result from the merging plasmoids rotating around each other.

9. 9.

   The plasmoid growth rate, regardless of growing in an isolated environment or through merging, is $\sim 0.1c$.

Our method clearly shows that plasmoids and mergers can be identified by multiple physical quantities (in our paper they are the magnetic vector potential and particle density). Although our segmentation is done manually, the identified plasmoids and mergers and their evolution can be fed into machine learning to segment plasmoids and mergers automatically. This can significantly reduce the computational cost and human labor involved in our present study, because right now we have to manually track the plasmoid evolution and merging processes case by case and code the envelops accordingly. In the future with machine learning, we can perform the above studies on more and longer combined PIC and radiation transfer simulations to boost the statistics of plasmoids and mergers, thus reducing the systematic uncertainties in modeling the radiation and polarization signatures from magnetic reconnection.

Our results are readily applicable to model 2D reconnection. However, three more steps are necessary to better study the multi-wavelength radiation signatures from magnetic reconnection. Firstly, we only study the synchrotron emission in this paper, using the optical band that is close to the cooling break as an example. Zhang et al. (2020) shows that synchrotron emission at different parts of the spectrum can behave differently in both the flux and polarization. Additionally, Zhang et al. (2022) finds that the synchrotron self Compton and external Compton exhibit different light curves compared to the synchrotron emission. In particular, the synchrotron self Compton can produce very fast flares due to the highly inhomogeneous particle and photon distributions in the secondary reconnection layer formed by plasmoid mergers. Although the present work can study all the above signatures, we do not have adequate statistics to reach any conclusive results. Therefore, we plan to investigate the multi-wavelength patterns with future machine learning studies. Secondly, our results are limited to 2D reconnection. 3D PIC simulations have shown that turbulence can become important in magnetic reconnection (Guo et al., 2021) and that substantial particle acceleration takes place very efficiently in the reconnection upstream (Giannios, 2010; Zhang et al., 2021b; Zhang et al., 2023). Although 2D slices along the guide field direction of the 3D simulations appear similar to 2D reconnection, the flux ropes, which are the 3D counterparts of plasmoids, can have much more complicated structures. Consequently, their mergers are much more complicated as well. Moreover, 3D simulations naturally allow different viewing angles, which can have a major impact on radiation and polarization. It is therefore necessary to examine that the radiation signatures obtained from 2D simulations are statistically consistent with those from the 3D simulations. Finally, PIC simulations cannot extend to realistic blazar zone scales with the available computing power, thus the initial and boundary conditions of the reconnection region may be very different from the realistic situation. Radiation signatures based on PIC simulations have to be tested against future hybrid simulations involving fluid dynamics, particle transport, and radiation transfer.