On the Rate of Crustal Failures in Young Magnetars

Neutron stars (NSs) are the final remnants of many massive stars after their supernova explosions. They manifest themselves thanks to the huge magnetic energy stored inside and around them. The detected emission can be powered by the electromagnetic torque that induces rotational energy losses, by accretion, by crustal heating, and by magnetic instabilities. In particular, magnetars, the most magnetic objects within the NS population, show magnetically driven enhancements of their thermal and nonthermal emission, called outbursts (Beloborodov & Li 2016; Coti Zelati et al. 2018). Their irregular and spectacular activity also includes short bursts and flares, believed to involve the magnetosphere, and possibly triggered by the interior dynamics.

Possibly related to magnetars (Margalit & Metzger 2018) are fast radio bursts (FRBs), very bright (∼50 mJy–100 Jy) millisecond-long bursts in the radio band, known for over a decade. Their distances as derived from the radio dispersion measures (DMs) have established their extragalactic nature, since the first identifications of their host galaxies confirmed their extragalactic nature (for recent reviews see, e.g., Cordes & Chatterjee 2019; Petroff et al. 2019). Their short durations, coherent radio emission, and the discovery of several repeating FRBs (Spitler et al. 2014; Fonseca et al. 2020) strengthened their connection with NSs, and in particular with strongly magnetized ones (Cheng et al. 2020; Wadiasingh et al. 2020).

Among the tens of theoretical models (see Platts et al. 2019 for a list), Lyutikov et al. (2016) propose that FRBs are giant pulses coming from young NSs spinning at very fast rotation frequencies. On the other hand, many works assume that the star's magnetic energy is the ultimate engine behind the emission. For instance, Lyubarsky (2014) proposes maser synchrotron emission arising from far regions, where the magnetar wind interacts with the nebula inflated by the wind within the surrounding medium, and generates a forward and a reverse shock, both able to provide intense, coherent emission frequencies of tens to hundreds of MHz. A variant is the baryonic shell model, originally proposed by Beloborodov (2017) and further developed by Metzger et al. (2019) and Margalit et al. (2020), in which an ultrarelativistic flare collides with matter ejected from an earlier flare. Kumar & Bošnjak (2020) propose a model where a relatively small Alfvén disturbance at the surface, with amplitude δB ∼ 10¹¹ G, is able to trigger the FRB emission via an efficient conversion to coherent curvature radiation from plasma bunches, accelerated by local electric fields.

The recent detection of an FRB-like radio burst simultaneous with a bright hard X-ray burst from the Galactic magnetar SGR 1935+2154 (Bochenek et al. 2020; Mereghetti et al. 2020; The CHIME/FRB Collaboration et al. 2020) showed for the first time that magnetar bursts can indeed be associated with FRB-like radio emission. It has been speculated that the underlying mechanism is the same as for the extragalactic events, for which we see only the very bright events from a much larger volume (Margalit et al. 2020).

Several works (e.g., Lu et al. 2020; Yuan et al. 2020) have proposed a magnetic disturbance from the surface driving the FRB. Regardless of the location and the exact mechanism, in the majority of magnetically powered FRB models, the ultimate trigger to the event is expected to come from the interior of the NS. Internal magnetic fields evolve in the long term due to the ohmic dissipation and Hall drift in the crust (Viganò et al. 2013; Pons & Viganò 2019). In the core, recent significant progress has been made in modeling and understanding the ambipolar diffusion and the evolving MHD equilibrium in superfluid (e.g., Graber et al. 2015; Bransgrove et al. 2017; Gusakov 2019) and non-superfluid cores (e.g., Castillo et al. 2017, 2020; Ofengeim & Gusakov 2018; Dommes et al. 2020), which for extreme values of magnetic field may proceed on relatively short timescales, thus changing the magnetic field on the crust–core boundary (Beloborodov & Li 2016).

In this Letter, we evaluate how often the crust of a magnetar will fail during the first centuries of its life, using an updated version of our magnetothermal code (Viganò et al. 2012). This study improves and extends to the early-stage previous results obtained for middle-aged magnetars (Perna & Pons 2011; Pons & Perna 2011), by including the Hall effect and up-to-date microphysics, but not the mechanisms in the core mentioned above. Since the early-stage dynamics are extremely sensitive to the initial conditions and the field topology, we perform the analysis under different assumptions on the initial magnetic field strength and geometry. We discuss our results in light of the current knowledge of FRB rates, recurrent times, and predictions in the FRB–magnetar models.

During the first decades of an NS's life, the star quickly cools from its initial 10¹¹–10¹² K to a few times 10⁸ K. As the temperature drops below the density-dependent melting value, the crust grows as a solid lattice of heavy ions. The crust freezing starts from the inner region a few minutes after birth, but takes several years to extend to the outer crust (e.g., Figure 8 of Aguilera et al. 2008). The crust is thought to provide an elastic response to stresses (Chugunov & Horowitz 2010) up to a maximum value

$$\begin{eqnarray}&&{\sigma }^{\max }=\left(0.0195-\displaystyle \frac{1.27}{{\rm{\Gamma }}-71}\right){n}_{i}\displaystyle \frac{{Z}^{2}{e}^{2}}{a},\end{eqnarray} \tag{ 1 }$$

where Γ = Z²e²/aT is the Coulomb coupling parameter, $a={\left[3/(4\pi {n}_{i})\right]}^{(1/3)}$ the ion sphere radius, n_i the ion number density, Z the charge number, e the elementary charge, and T the temperature. A recent calculation (Kozhberov & Yakovlev 2020) confirms the validity of the approximation (Equation (1)).

As the magnetic field evolves, the magnetic tensor ${M}_{{ij}}({\boldsymbol{x}},t)={B}_{i}({\boldsymbol{x}},t){B}_{j}({\boldsymbol{x}},t)/(4\pi )$ changes. Moreover, a stress quantified as a difference between M_ij and an equilibrium ${M}_{{ij}}^{\mathrm{eq}}$ piles up in the crust. When it reaches the maximum stress locally allowed, $({M}_{{ij}}-{M}_{{ij}}^{\mathrm{eq}})\simeq {\sigma }^{\max }({\boldsymbol{x}})$, the crust fails. After the event, a new equilibrium state is established, the affected region freezes again, and the crust starts to respond elastically, storing further stress until the cycle is repeated.

An evaluation of the frequency of such events, which can potentially trigger magnetar activity, can be done by properly following the magnetothermal evolution and the accumulated local stresses. Such simulations have been performed by Perna & Pons (2011) and Pons & Perna (2011), considering the age range of the observed sources, ∼400–10⁵ yr. They used the magnetothermal code by Pons et al. (2009), not considering Hall effects.

As in Pons & Perna (2011) and Perna & Pons (2011), we assume that after a crustal failure event, the surrounding regions of the crust, being close to the maximum stress, M_ij ≥ σ^max, are affected by the instability and resettle to equilibrium. In our simulations, we set the effective parameter to  = 0.9. Note that increasing it to values closer to 1 will lead to more frequent events but smaller involved regions, and vice versa.

Our simulations present a few relevant simplifications: (i) we do not include self-consistent magnetoelastic evolution of the crust (Cumming et al. 2004; Li et al. 2016; Bransgrove et al. 2017; Thompson et al. 2017), and the local magnetic field reconfiguration likely produced after a crustal failure; (ii) in the core we only consider the ohmic term, which acts on timescales much longer than in the crust: in practice, the magnetic field does not evolve in the kyr here considered and a current sheet develops at the core surface; (iii) at the NS surface we impose a potential magnetic field solution (as in almost all studies, reviewed by Pons & Viganò 2019).

The main challenge in considering early ages is that the results are extremely sensitive to the initial magnetic field configuration, for both physical and numerical reasons. Before freezing, the magnetized proto-neutron star fluid will settle to an MHD equilibrium/stationary state. Several works have proposed initial configurations, most of them based on a very large-scale and smooth initial magnetic field, usually with a purely dipolar poloidal component. However, nature typically provides much more turbulent, off-axis and nonsymmetric configurations, as planetary and stellar magnetic fields show. Implementing such complicated configurations in simulations is a hard task. First of all, the possible solutions are potentially infinite and it is not clear which ones could be more realistic. Second, having small scales means having much faster Hall drift and whistler waves, characteristic of the Hall MHD equations that govern the solid crust. The propagation of such relatively fast transient waves dominates the dynamics and the piling up of the stresses over the first centuries. However, simulating such small-scale waves requires very fine space and time resolution, making the computational cost unfeasible, and the numerical instabilities become hardly controllable.

3.1. Initial Configurations

Aiming to reduce numerical instabilities that could provoke artificial crustal failures, we have used the latest version of the finite-volume axially symmetric magnetothermal code, improved in stability and efficiency compared to Viganò et al. (2012, 2013). It includes numerical methods that are suitable to the finite-volume integral version of the Hall induction equation, and are optimized in terms of stability, accuracy, and efficiency. We also employ the most updated microphysical inputs (for a review of cooling and transport see Potekhin et al. 2015). We make sure that results are close to the numerical convergence (we have variations of the event rate up to a maximum of 30% for different resolutions, much less than the intra-model variability), using a spatial resolution of 100 × 200 points in the angular and radial directions, and a fixed time step of ∼10⁻⁴–10⁻³ yr for the magnetic advance, and of 10⁻² yr for the temperature evolution.

To assess the sensitivity of results upon the uncertain initial conditions, we have considered very different topologies. Although not necessarily realistic, the variety considered allows us to explore the range of crustal failure events. They are summarized in Table 1, together with the evolution of the events' frequency during the first 1 kyr of the NS's life.

Table 1.  Initial Models Considered and Events Simulated

Model	Topology	Bdip	${E}_{\mathrm{mag}}^{\mathrm{cr}}$	${E}_{\mathrm{tor}}^{\mathrm{cr}}/{E}_{\mathrm{mag}}^{\mathrm{cr}}$	N100	N400	N1000	First Event
		(G)	(erg)					(yr)
CrDip	CC: PD+TQ	1014	3.2 × 1046	∼50%	107	192	245	1.3
CrDipL	CC: PD+TQ	1013	3.2 × 1044	∼50%	0	0	0	⋯
CrDipH	CC: PD+TQ	3 × 1014	2.9 × 1047	∼50%	2380	3315	3724	0.4
CrDipE	CC: PD+TQ	1015	3.2 × 1048	∼50%	23360	30842	34890	0.1
CrMultiL*	CC: PM+TM	1013	9.3 × 1044	∼50%	28	41	49	3.0
CrMultiH*	CC: PM+TM	1014	2.9 × 1047	∼50%	15538	21397	25513	0.1
CrMultiT*	CC: PM+TM	1013	4.8 × 1046	∼99%	525	1335	1908	0.7
CrDipD	CC: PD+TD	1014	3.4 × 1046	∼50%	88	176	254	1.6
CrDipP	CC: PD+TQ	1014	1.6 × 1046	∼1%	45	58	67	2.6
CoDipL	CT: PD+TD	1014	4.7 × 1044	∼50%	10	10	10	13.7
CoDipH*	CT: PD+TD	3 × 1014	8.8 × 1045	∼50%	160	168	174	1.6
CoMulti	CT: PM+TD	1014	7.5 × 1046	∼50%	1948	2265	2370	0.8
CoMultiP	CT: PM+TD	1014	2.6 × 1046	∼1%	528	577	597	4.0
CoMultiPH*	CT: PM+TD	3 × 1014	2.4 × 1047	∼1%	9497	9997	10460	0.8
CoMultiPE*	CT: PM+TD	6 × 1014	9.4 × 1047	∼1%	21090	23840	29065	0.4

Note. Topology (where CT is a core-threaded configuration, CC is a crust-confined configuration, P and T are the poloidal and toroidal field, D and Q mean purely dipolar and quadrupolar, and M is a mixed multipolar configuration); surface value of the dipolar poloidal component B_dip; magnetic energy ${E}_{{\rm{mag}}}^{{\rm{cr}}}$ stored in the crust, the fraction ${E}_{{\rm{tor}}}^{{\rm{cr}}}/{E}_{{\rm{mag}}}^{{\rm{cr}}}$ stored in the toroidal field; N^x is the number of events shown during the first x = 100, 400, and 1000 yr; first event is the age at which the crust fails for the first time. Models marked with * show small local numerical instabilities that can provide artificial events, quantifiable by up to a maximum of 30% of the shown values N, thus not affecting the order of magnitude.

Download table as:  ASCIITypeset image

On one hand, we take into account crust-confined models with different weights in the initial multipoles. While it is unlikely that the field is completely expelled from the core, such a kind of topology can be useful to provide an upper limit on the frequency of outburst/burst activity triggered by the dynamics, since the field lines are not tied to the highly conductive core and the currents supporting the magnetic field flow entirely in the crust, where timescales are reduced compared to the core. Within these models, we vary the relative fraction of magnetic energy stored in different multipoles of both the poloidal and toroidal components, as follows:

• 1.  
  In the fiducial model CrDip we consider equipartition of magnetic energy between a dipolar poloidal magnetic field with polar surface value B_dip = 10¹⁴ G and a toroidal quadrupolar field having the same magnetic energy as the poloidal one (and a maximum value ${B}_{\varphi }^{\max }=1.27\times {10}^{15}$ G).
• 2.  
  In the low, high, and extreme field models CrDipL, CrDipH, and CrDipE, we assume the same topology but with the magnetic field intensity rescaled by a factor 0.1, 3, and 10, respectively.
• 3.  
  In the low and high field models CrMultiL and CrMultiH, we consider a topology with a similar contribution of the first four multipoles in both the poloidal and toroidal components, with an equipartition of magnetic energy between them. In CrMultiT, the quadrupolar toroidal field accounts for 99% of the total magnetic energy.
• 4.  
  In model CrDipD a dipolar (instead of a quadrupolar) toroidal field is considered.
• 5.  
  In model CrDipP we adopted a poloidal-dominated configuration, with the quadrupolar toroidal field accounting for less than 1% of the total magnetic energy.

Additionally, we consider different core-threaded configurations:

• 1.  
  Model CoDipL (low field) and CoDipH (high field), consisting of the often used large-scale twisted dipole plus torus configurations, as in Akgün et al. (2017). Again, we set an equal fraction of poloidal and toroidal magnetic energy in the crust.
• 2.  
  In model CoMulti, we consider a topology with similar contributions from the dipole and the quadrupole in the poloidal field, and a dipolar toroidal field, with an equipartition of magnetic energy between them. On the other hand, for CoMultiP, CoMultiPH, and CoMultiPE, the toroidal magnetic energy accounts for less than 1% of the total magnetic energy, where H and E refer to high and extreme field models, respectively.

Note how models CrDip, CrMultiH, CrDipD, CrDipP, CoDipL, CoMulti, and CoMultiP have the same dipolar polar magnetic field B_dip = 10¹⁴ G, which rules the spin-down and, indirectly, the rotationally powered emission. However the magnetic energy stored in the crust, ${E}_{\mathrm{mag}}^{\mathrm{cr}}$, can vary by a few orders of magnitude, depending on the relative weight of the multipolar and toroidal components and on whether currents circulate only in the crust or not.

In Figure 1 we show the initial and evolved (at 1 kyr) magnetothermal models CrDip (top panels), CrMultiH (central panels), and CoDipL (bottom panels). The crustal-confined models show relevant magnetic activity (left hemispheres of each panel), with a local nontrivial rearrangement of the toroidal field distribution (colors) and a deformation of the poloidal field lines (white lines). The temperature distribution (right hemisphere) also evolves to be inhomogeneous, due to the transport anisotropy induced by the local intensity and direction of magnetic field. On the other hand, the core-threaded model barely shows any evolution, due to the fact that the currents are mostly concentrated in the core, where the conductivity is very high and no Hall effect is present: the magnetic field lines are rooted and almost frozen in the core. The core-threaded cases give much slower dynamics than the crust-confined ones, due to the large characteristic spatial scales, and the fact that most of the supporting currents circulate in the highly conductive, non-solid core. A realistic configuration could be core-threaded, but consists of small scales, both in the core and the crust, so that the corresponding results should be effectively covered by our exploration of parameters. Note that, in all cases, the value of the dipolar component at the surface, B_dip, decreases at most by a few percent.

Zoom In Zoom Out Reset image size

3.2. Event Rate

Focusing on the predicted event rate, the gross number for model CrDip is ∼250 events during the first 1 kyr of the NS's life, mostly concentrated in the first century. Furthermore, we emphasize that by changing the overall magnetic strength with respect to model CrDip, the frequency decreases for the low field dipolar model (CrDipL), and increases for the high field (CrDipE) and the extreme field (CrDipH) model, by ∼1 and ∼2 orders of magnitude, scaling the number of events as ${N}^{1000}\sim {E}_{\mathrm{mag}}^{\mathrm{cr}}/{10}^{44}$. The same scaling is noticed for the core-threaded models. There are no significant differences between a different topology for the toroidal field (quadrupolar in CrDip and dipolar in CrDipD), while important differences are noted if the magnetic energy is changed, maintaining the same B_dip: poloidal-dominated CrDipP, or considering multipolar configurations (CrMultiL, CrMultiH, and CrMultiT).

For all models, the older the magnetar, the fewer the expected number of events. As shown in Table 1, ∼40%–85% of the events during the first 1000 yr are concentrated in the first century, with a gradual decrease. This percentage decreases to ∼10%–40% in the next 300 yr, whereas in the last 600 yr, the percentage oscillates between ∼5% and 20%. This is consistent with the findings for middle-aged magnetars (Perna & Pons 2011; Pons & Perna 2011): magnetars become less and less active, but they can still show activity even at late stages.

The waiting time distributions are illustrated in Figure 2. In the left panel, we compare models CrDip, CrDipL, CrDipH, and CrDipE. We find that increasing the strength of the magnetic field results in a shorter waiting time. This is also the case if we consider a multipolar topology (models CrMultiH and CoMulti), whereas lowering the fraction of the toroidal energy (model CrDipD) has a tiny impact on the distribution. However, the waiting time is shifted from ∼[10:10⁴] to ∼[0.1:300] days, by increasing the field strength and considering the multipolar topology.

Zoom In Zoom Out Reset image size

In Figure 3, we present the crustal failure distributions as a function of the NS radius for a set of models, i.e., models CrDip, CrDipE, and CoMulti. Even though such distributions are partially model-dependent, most of the events happen in the outer crust, i.e., R ∼ [11.3:11.55] km, because the maximum stress is smaller in the outermost, lighter layers. A nonnegligible fraction of the events takes place in the outermost 50 m of the NS, i.e., ∼25%.

Zoom In Zoom Out Reset image size

Ideally, modeling the link between magnetospheric magnetar activity and crustal failure requires (i) following the internal magnetic dynamics and how the magnetic stresses grow with time, (ii) understanding and quantifying the mechanisms leading to the propagation of the disturbance into the magnetosphere, and (iii) describing how the magnetospheric instability/reconnection powers the emission.

This study focuses on the first issue, assessing the crustal failure frequency of newly born magnetars, under different initial topologies, through 2D magnetothermal simulations, including the Hall term. Moreover, the presence of small-scale structures, very likely if the birth amplification involves turbulent processes, can increase the expected event rate, through initial enhanced Hall dynamics. More factors can further extend the range of events found, among which we must mention the magnetic field boundary conditions (considered here as potential, possibly inadequate for the dense environment around a newborn magnetar) and the fixed effective parameter , controlling how large the region is undergoing failure at the same time.

The main finding is that the event frequency scales linearly with the magnetic energy stored in the crust, with some intrinsic dispersion given by the topology, ${E}_{\mathrm{mag}}^{\mathrm{cr}}$ (see Figure 4), which strongly depends on the assumed and highly uncertain field topology.

Zoom In Zoom Out Reset image size

Given the simplifications in our approach, additional words of caution are needed in the interpretation of the results:

(i) We do not include the recent advances in the modeling of the fluid motion in the NS core (with or without superfluidity), limiting to either confine the field to the crust, or to consider only the very slow ohmic dissipation. In particular, the ambipolar diffusion could drive (slightly) faster changes and higher rates (Beloborodov & Li 2016).

(ii) Magnetoelastic simulations, and the inclusion of the local field rearrangement and dissipation, could reduce the rates, especially the most extreme ones (thus potentially weakening the linear trend in the upper range). Moreover, the critically stressed crust may behave plastically (Kobyakov & Pethick 2014), leading to Hall drift and thermoplastic waves (Beloborodov & Levin 2014; Li et al. 2016) and partially hampering the Hall dynamics (Lander & Gourgouliatos 2019).

(iii) 3D dynamics and more realistic boundary conditions could also have a quantitative impact, difficult to guess a priori.

We expect that overcoming such caveats above will systematically change the quantitative results, but not the important conclusion that the dipolar field strength is a barely relevant parameter when considering a young magnetar's activity. As a matter of fact, for a given B_dip, our calculations show a very broad range of crustal event rates. Even a relatively modest B_dip = 10¹⁴ G could potentially provide almost no events (large-scale field penetrating in the core) or N ∼ 1 event per day, if most of the magnetic energy is hidden in the form of crustal multipolar and toroidal fields. The expected crustal failure rate is, thus, strongly dependent on the initial topology of the magnetic field. These results shed light on a potential issue of most FRB–magnetar models, which are usually focused on the value of B_dip only.

4.1. Comparison with Fast Radio Bursts and Galactic Magnetars

Our calculations have implications both for FRBs (if indeed associated with magnetars), as well as for Galactic magnetars. In both scenarios, the observed event rate can be summarized as

$$\begin{eqnarray}&&{N}_{\mathrm{obs}}={N}_{\mathrm{failure}}\times {N}_{\mathrm{sources}}\times {x}_{\mathrm{vis}},\end{eqnarray} \tag{ 2 }$$

where N_sources is the number of potential sources (young magnetars within a certain observable volume), N_failure is the crustal event rate per star calculated in our study, and x_vis is the effective fraction of events that will actually be detected. The latter number is basically unconstrained a priori, depending on the specific mechanisms ultimately leading to a detection, to observational biases, and to both physical (beaming, distance, absorption, etc.) and instrumental (sensitivity, field of view, bandwidth, etc.) limits.

In the case of magnetar bursting activity happening in the magnetosphere, x_vis depends on the crustal response to failure (one can expect that only the outermost events can further affect the magnetosphere), the propagation of the disturbance outside, the magnetospheric dynamics and the emission mechanism. Besides these fundamental physical uncertainties, emission beaming, distance, and absorption will affect x_vis, resulting in an expected value much smaller than 1.

For FRBs, according to Rane et al. (2016), the observed rate is estimated to be N_obs ∼ 10³–10⁴ per day per sky above 4 mJy per ms. Considering as a reference the number of galaxies within the distance of the prototypical FRB repeater FRB 121112, ∼10⁷, and assuming ∼10 magnetars per galaxy with ≲1000 yr, we would expect a potential number of FRB-emitting magnetars of about N_sources ∼ 10⁸. Equating these numbers, we obtain a constraint ${x}_{\mathrm{vis}}^{\mathrm{frb}}\times {N}_{\mathrm{failure}}\sim {10}^{-5}\mbox{--}{10}^{-4}$ per day. This means that models with high crustal magnetic energy, ${E}_{\mathrm{mag}}^{\mathrm{cr}}\gtrsim {10}^{46}$ erg, which have N_failure ∼ 1 per day, are compatible with x_vis ≪ 1.

Similarly, very young (≲1 kyr) Galactic magnetars are expected to be N_sources ∼ 10. The observed rate of outbursts for them (excluding older magnetars) is N_obs ∼ 10⁻³–10⁻⁴ per day (Coti Zelati et al. 2018). Therefore, ${x}_{\mathrm{vis}}^{\mathrm{galactic}}\times {N}_{\mathrm{failure}}\sim {10}^{-5}\mbox{--}{10}^{-4}$, similarly to the FRBs.

This preliminary estimate, although plagued by uncertainties on the total number of sources and the still incomplete statistics, shows that the range of events detected per star, x_vis × N_failure, in the FRB and Galactic magnetar outburst scenarios, is at least compatible. This points to a possible common origin of crustal triggers, even though they manifest in different ways.

Our simulations show that the dipolar component, on which most FRB–magnetar models are currently based, plays a minor role in determining the bursting activity of the magnetosphere if triggered by crustal failures. Providing reliable numbers for the expected event rate is intrinsically difficult, mainly for the intrinsic limitations of the approach and the highly uncertain initial configuration of the magnetic field. Here we have shown how, by varying the configuration from a large-scale, smooth core-threaded field, to a multipolar field confined in the crust, the calculated crustal failure rate changes by several orders of magnitude, with the same dipolar poloidal value, B_dip. However, we have found that the crustal magnetic energy is instead a good tracer of the expected number of events, regardless of the exact topology.

Further observations and a better understanding of the mechanisms leading to the propagation of the disturbance and the emission of bursts will help constrain the fraction of events that can be detectable (x_vis). If this parameter was constrained even at an order-of-magnitude level, it would allow us to constrain the initial topology of the magnetic field.

We thank Alice Borghese, Francesco Coti Zelati, and Vanessa Graber for helpful suggestions. C.D., D.V., N.R., and A.G.G. are supported by the ERC Consolidator Grant "MAGNESIA" (No. 817661) and acknowledge funding from grants SGR2017-1383 and PGC2018-095512-BI00. J.A.P. acknowledges support by the Generalitat Valenciana (PROMETEO/2019/071) and by AEI grant PGC2018-095984-B-I00. R.P. acknowledges support from NSF award AST-1616157. We acknowledge support from the PHAROS COST Action (CA16214).