© The Authors 2025

1 Introduction

Similar to external galaxies, the Milky Way is permeated by magnetic fields. Magnetic fields are relevant in various processes and scales, from regulating star formation to shaping large-scale structures in the disk and the halo of the Galaxy. Magnetic fields in the Galaxy also affect our ability to study the Universe’s structure and history and to identify the sources of ultra-high energy cosmic rays (e.g., see Beck 2015; Haverkorn 2015 for recent reviews).

To gain insight into the general structure and amplitude of the Galactic magnetic field (GMF) on the kiloparsec scale, authors have primarily relied on observations of Faraday rotation measures (RM) toward extragalactic sources, Galactic synchrotron emission, and polarized emission from Galactic dust. Using forward modeling frameworks, where model predictions are compared to observations, authors have deduced constraints on a variety of three-dimensional (3D) models for the GMF (e.g., see Jaffe 2019 for a recent review on model tuning). One of the main obstacles to such an approach is the line-of-sight-integrated nature of the observables, which leads to degeneracy between models and model parameters and, consequently, uncertainties in the reconstructed picture of the GMF. In the future, this degeneracy will probably be overcome, at least for the first few kiloparsecs around the Sun, by the advent of tomography data derived from both starlight polarization and Faraday tomography (e.g., Pelgrims et al. 2023, 2024; Erceg et al. 2024a,b), which are made possible by major surveys (Magalhães et al. 2005; Tassis et al. 2018; Shimwell et al. 2022). In the meantime, it is difficult to disentangle the respective contributions of the Galaxy’s near and far regions to the line-of-sight-integrated observables.

In this context, it should be stressed that most efforts to model the large-scale GMF have overlooked the fact that the Solar System resides in the Local Bubble. The Local Bubble is a cavity with an unusually low density of gas (e.g., Cox & Reynolds 1987; Lallement et al. 2003), filled with an X-ray emitting plasma (e.g., Liu et al. 2017), and surrounded by a shell of cold (ionized) gas and dust (e.g., Pelgrims et al. 2020; O’Neill et al. 2024b). This local cavity has most likely been created by successive supernova explosions that occurred in the past 10– 15 Myr (e.g., Maíz-Apellániz 2001; Breitschwerdt et al. 2016; Schulreich et al. 2023), and currently has an irregular shape, with a radius of around 100 pc in the Galactic plane and around 300 pc perpendicular to it. As for other bubbles in the interstellar medium (ISM), it is expected that supernova explosions and the resulting motion of the swept-up matter have deformed the initial magnetic field and amplified its strength in compressed regions. This is in agreement with the picture obtained from starlight polarization that the magnetic field in the local ISM does not follow the large-scale GMF (e.g., Heiles 1998; Leroy 1999; Santos et al. 2011; Frisch et al. 2012; Berdyugin et al. 2014; Gontcharov & Mosenkov 2019) and that the amplitude of the magnetic field is likely enhanced in the shell of the Local Bubble (Andersson & Potter 2005, 2006; Medan & Andersson 2019).

Although imprints of the magnetized shell of the Local Bubble have recently been detected in tomographic data, from measurements on Faraday depth and starlight polarization, in combination or not with dust polarized emission (Erceg et al. 2024a,b; Skalidis & Pelgrims 2019; Pelgrims et al. 2024; Uppal et al. 2024), its contribution to line-of-sight-integrated observables which are sensitive to the GMF is still unknown, or at least not very constrained. The chance for its contribution to be dominant naturally increases toward high-Galactic latitudes as less matter is encountered beyond the Local Bubble wall. This is even more likely to be the case the smaller the scale height of the density distribution of the relevant ISM matter. This, for example, motivated Alves et al. (2018) to use the dust polarized emission at high-Galactic latitudes to constrain a physical model for the magnetic field in the shell of the Local Bubble whose shape was assumed to be spheroidal, and with no thickness. This motivation was later validated by Skalidis & Pelgrims (2019) (but see Halal et al. 2024) and encouraged Pelgrims et al. (2020) to rework the analysis of Alves et al. (2018) with a more realistic shape for the Local Bubble shell that they deduced directly from 3D dust density maps available at the time. However, it is difficult to predict how the relative contributions from the Local Bubble shell and from the background change with sky regions, toward lower Galactic latitudes. This is primarily due to depolarization effects that occur both along the sightline and within the beam of observation. Surprisingly, perhaps, O’Neill et al. (2024a) suggested recently that the wall of the Local Bubble may be a substantial contributor to the dust-polarized emission for large portions of the sky.

Recently, Korochkin et al. (2025) were the first to include a prescription for the magnetic field in the shell of the Local Bubble into the modeling of the regular component of the GMF using full-sky data on Galactic synchrotron polarized emission and Faraday RM from extragalactic sources. They found that the contribution from the shell of the Local Bubble is substantial. However, in their implementation, they used a spherical model for the shape of the Local Bubble and a geometrical prescription to derive the magnetic field in its shell, which did not lead to a divergence-free magnetic field. With this paper, we aim to go beyond this and quantify the effects of a more realistic shape for the shell of the Local Bubble.

In their paper, Alves et al. (2018) derived an analytical expression for a divergence-free magnetic field in the shell of the Local Bubble that they assumed to be very thin. While in principle this formalism makes it possible to describe the magnetic field in non-spherical shells, it has the important drawback that the shell must be infinitely thin and, as a consequence, that the amplitude of the magnetic field becomes very large (infinite) in the shell. Consequently, this analytical prescription is inadequate for characterizing the magnetic field in a realistic shell of finite thickness, and for predicting the shell’s contribution to Faraday RM and synchrotron, both of which depend on the amplitude of the magnetic field. In this work, we move past the assumption that the bubble shells must be infinitely thin. We propose a general analytical expression for the divergence-free magnetic field in the thick shell of bubbles which, additionally might be nonspherical. The paper is structured as follows.

We derive our analytical model in Sect. 2 where we also present its main features and some illustrative examples. In Sect. 3, we particularize our model to the case of the Local Bubble and estimate its contribution to the full-sky data on synchrotron polarized emission and Faraday RM, adopting realistic shapes for the shell of the Local Bubble as inferred from a 3D dust map. We discuss the limitations of our model and the possible implications of our findings for the modeling of the large-scale GMF in Sect. 4.

2 An analytical model for the magnetic field in thick bubble shells

In this paper, we build upon the pioneering work of Parker (1970) and Alves et al. (2018) to describe the divergence-free magnetic field in thick shells of ISM bubbles. We only consider the coherent (regular) component of the magnetic field. We first derive the equations in a general case without specifying any functional form for the displacement vector field that relates the presenttime position of test particles to their initial locations. Then, we particularize our solution to the case where the displacement vector field is linear in the radial coordinate (with respect to the center of explosion) and explore our model predictions for various possible geometries for the bubble shell.

2.1 General expression for the magnetic field

To derive the equations for the magnetic field in the thick shell of a bubble, we closely follow the derivation in Alves et al. (2018) that is based on the work of Parker (1970). Our base assumptions are that an ISM bubble results from a single explosion that has swept up matter into a shell with finite thickness and that the explosion-induced motions are purely radial with respect to the explosion center, but do not necessarily have a spherical symmetry. That is, the explosion-induced velocity flow is not necessarily isotropic. We further adopt the frozen-in approximation that implies that magnetic field lines follow the matter in its motion. As a result of this coupling, the vector field of the magnetic field, which we assume to be initially uniform in strength and direction, is modified by the explosion. The magnetic field lines are deformed to closely follow the geometry of the shell and are confined within it, thus increasing the amplitude of the magnetic field. The amplitude of the magnetic field vanishes inside the bubble and remains unchanged outside it.

To proceed, we consider that the explosion center is at the origin of a Cartesian coordinate system (x, y, z) to which is associated a spherical coordinate system (r, θ, ϕ) with its corresponding unit basis vectors (e_r, e_θ, e_ϕ). In this reference frame, the velocity field at any time t takes the form V = V(r,θ,ϕ) e_r, and the initial position of a test particle presently at position r = r e_r can be written as $$\begin{array}{l}{r}^0=r^0(r,\theta ,\varphi )e_r\\ =r-\Lambda ,\end{array}$$(1)

where we introduced the displacement vector field Λ that results from the time integration of V. As it is the case for V, the displacement vector field is radial but not necessarily with spherical symmetry. We have Λ = Λ(r, θ, ϕ) e_r.

Generally, the frozen-in approximation is given by the induction equation in the limit of zero magnetic diffusion. Therefore, the magnetic field B must be given by the following: $$\frac{\partial B}{\partial t}=\nabla \times (V\times B).$$(2)

Introducing the vector potential A, such that B = ∇ × A, in Eq. (2), swapping the time and spatial derivatives, and decurling, we obtain $$\frac{\partial A}{\partial t}=V\times (\nabla \times A)+\nabla f,$$(3)

where f is an arbitrary scalar field. Adopting the gauge f = −V · A, we then write the equation for the Cartesian components of A as follows: $$\begin{array}{l}\left[\frac{\partial }{\partial t}+{\displaystyle \sum _{j=1}^3{V}_j}\frac{\partial }{\partial x_j}\right]A_i=-{\displaystyle \sum _{j=1}^3{A}_j}\frac{\partial V_j}{\partial x_i}\\ \iff \frac{ \text{d}}{ \text{d}t}{A}_i=-{\displaystyle \sum _{j=1}^3{A}_j}\frac{\partial V_j}{\partial x_i},\end{array}$$(4)

where the operator on the left-hand side is known as the advective derivative. Integrating in time the above equation leads to the expression of the Cartesian components of the potential vector in the present time and at location (x) as a function of the initial potential vector at the initial location (x⁰). We have $$A_i(x)={\displaystyle \sum _{j=1}^3{A}^0}j\left(x^0\right)\frac{\partial x_j^0}{\partial x_i}.$$(5)

This is Eq. (10) in Parker (1970). IfA⁰(x⁰) is known, determining A(x) amounts to determining how the moving coordinate system evolves in space with respect to the fix one. Since we assume a radial velocity field, the displacement vector field is also radial, and we can write its Cartesian components as $${(}_{\Lambda }=\text{Λ}\frac{\partial x_i}{\partial r},$$(6)

and therefore $$x_i^0=x_i-\text{Λ}\frac{\partial x_i}{\partial r}.$$(7)

Inserting this into Eq. (5) leads to $$A_i={{\displaystyle \sum _{j=1}^3{A}^0}}_j{\delta }_{ij}\left(1-\frac{\text{Λ}}{r}\right)-{{\displaystyle \sum _{j=1}^3{A}^0}}_j\frac{\partial x_j}{\partial r}\left(\frac{\partial \text{Λ}}{\partial x_i}-\frac{\text{Λ}}{r}\frac{\partial x_i}{\partial r}\right),$$(8)

where δ_ij is the Kronecker symbol. Then, using the fact that the spherical component of A can be written as $$\begin{array}{l}{A}_r={\displaystyle \sum _{i=1}^3{A}_i}\frac{\partial x_i}{\partial r}={\displaystyle \sum _{i=1}^3{A}_i}\frac{x_i}{r}\\ A_{\theta }={\displaystyle \sum _{i=1}^3{A}_i}\frac{1}{r}\frac{\partial x_i}{\partial \theta }\\ A_{\varphi }={\displaystyle \sum _{i=1}^3{A}_i}\frac{1}{r\sin \theta }\frac{\partial x_i}{\partial \varphi },\end{array}$$(9)

and that the same holds true for A⁰ (using the same spherical coordinates), we finally write the spherical components of A(r) as a function of A⁰(r⁰) and the displacement vector and its partial derivatives: $$\begin{array}{l}{A}_r=A_r^0\frac{\partial r^0}{\partial r}=A_r^0\left(1-\frac{\partial \text{Λ}}{\partial r}\right)\\ A_{\theta }=A_{\theta }^0\frac{r^0}{r}-\frac{1}{r}\frac{\partial \text{Λ}}{\partial \theta }{A}_r^0\\ A_{\varphi }=A_{\varphi }^0\frac{r^0}{r}-\frac{1}{r\sin \theta }\frac{\partial \text{Λ}}{\partial \varphi }{A}_r^0.\end{array}$$(10)

These equations are equivalent to Eqs. (A4) and (A5) derived in (Alves et al. 2018).

Now we can derive the spherical-coordinate components of B(r) simply by taking the curl of A(r). In spherical coordinates, we have $$\begin{array}{l}{B}_r=\frac{1}{r\sin \theta }\left(\frac{\partial \left(\sin \theta A_{\varphi }\right)}{\partial \theta }-\frac{\partial A_{\theta }}{\partial \varphi }\right)\\ B_{\theta }=\frac{1}{r}\left(\frac{1}{\sin \theta }\frac{\partial A_r}{\partial \varphi }-\frac{\partial \left(r{A}_{\varphi }\right)}{\partial r}\right)\\ B_{\varphi }=\frac{1}{r}\left(\frac{\partial \left(r{A}_{\theta }\right)}{\partial r}-\frac{\partial A_r}{\partial \theta }\right).\end{array}$$(11)

The same holds true for B⁰(x⁰) in the vector basis associated to the moving coordinate system (r⁰, θ⁰, ϕ⁰). Because the displacement vector field is a function of (r, θ, ϕ), the change of coordinates r ↔ r⁰ (or equivalently x ↔ x⁰), implies the following relations between the partial derivatives: $$\begin{array}{l}\frac{\partial }{\partial r}=\left(1-\frac{\partial \text{Λ}}{\partial r}\right)\frac{\partial }{\partial r^0}\\ \frac{\partial }{\partial \theta }=\frac{\partial }{\partial {\theta }^0}-\frac{\partial \text{Λ}}{\partial \theta }\frac{\partial }{\partial r^0}\\ \frac{\partial }{\partial \varphi }=\frac{\partial }{\partial {\varphi }^0}-\frac{\partial \text{Λ}}{\partial \varphi }\frac{\partial }{\partial r^0}.\end{array}$$(12)

The fact that Λ is radial, however, guarantees that the angular coordinates of a test particle remain unchanged: θ = θ⁰ and ϕ = ϕ⁰. Taking this into account, the analytical equations for the components of the magnetic field in the shell of a bubble become: $$\begin{array}{l}{B}_r={\left(\frac{r^0}{r}\right)}^2{B}_r^0+\frac{r^0}{r}{\nabla }_t\text{Λ}\cdot B_t^0\\ B_{\theta }=\frac{r^0}{r}\frac{\partial r^0}{\partial r}{B}_{\theta }^0\\ B_{\varphi }=\frac{r^0}{r}\frac{\partial r^0}{\partial r}{B}_{\varphi }^0,\end{array}$$(13)

where ∇_tΛ and $B_t^0$ are the orthoradial components of the gradient of Λ(r, θ, ϕ) and B⁰, that is, $${\nabla }_t=\frac{1}{r}\frac{\partial }{\partial \theta }{e}_{\theta }+\frac{1}{r\sin \theta }\frac{\partial }{\partial \varphi }{e}_{\varphi },\text{ and}$$(14) $$B_t^0=B_{\theta }^0{e}_{\theta }+B_{\varphi }^0{e}_{\varphi }.$$(15)

The expressions in Eq. (13) are general and have been obtained without the need to specify any form for the displacement vector field except that we required it to be radial with respect to the explosion center. The present-time magnetic field in the shell of the bubble can be obtained as soon as the initial magnetic field, the displacement vector field, and its partial derivatives are known. The expressions in Eq. (13) show that the magnetic field lines remain in the same plane defined by the explosion center and the orthoradial component of the initial magnetic field $B_t^0$, as measured in spherical coordinates centered on the explosion center. The radial component (B_r) adjusts so that magnetic field lines follow the bubble shape. It is worth noticing that, in general, the radial component of the magnetic field does not vanish in the bubble shell. The main challenge in estimating Eq. (13) resides in determining the partial derivatives of Λ(r, θ, ϕ) at any location within the bubble shell.

2.2 A linear model for the displacement vector field

To move forward, we need to particularize our equations with an expression for the displacement vector field Λ(r, θ, ϕ), that could ideally be constrained from observation. We consider that Λ(r,θ,ϕ) is linear in r and is such that all the matter located within the bubble at the initial time has been swept up in a shell defined by its inner and outer radius surfaces R_min and R_max , respectively. By construction, we assume that R_min and R_max are both functions of the angular coordinates (θ, ϕ) measured in the explosion-centric reference frame, and that the volume encompassed by the bubble is determined by r ≤ R_max. With these choices, the expression for the displacement along any radial direction may be written as: $$\text{Λ}(r,\theta ,\varphi )=\frac{R_{\min }}{R_{\max }-R_{\min }}\left(R_{\max }-r\right),$$(16)

where the angular dependence of R_min and R_max are implicit. In this simple case, the partial derivatives of Λ become: $$\frac{\partial \text{Λ}}{\partial r}=-\frac{R_{\min }}{R_{\max }-R_{\min }}$$(17) $$\begin{array}{l}\frac{\partial \text{Λ}}{\partial \theta }=r\frac{R_{\min }{\partial }_{\theta }{R}_{\max }-R_{\max }{\partial }_{\theta }{R}_{\min }}{{\left(R_{\max }-R_{\min }\right)}^2}\\ +\frac{R_{{\max }^2}{\partial }_{\theta }{R}_{\min }-R_{{\min }^2}{\partial }_{\theta }{R}_{\max }}{{\left(R_{\max }-R_{\min }\right)}^2},\end{array}$$(18)

where we have used the shorthand notation for the partial derivative: ∂_θ ≡ ∂/∂θ. ∂Λ/∂ϕ is given by an equivalent equation, replacing θ by ϕ in Eq. (18). As such, the partial derivatives of Λ in the shell of the bubble are fully determined by the shape of the inner and outer surfaces of the bubble shell (R_min , R_max , and their angular derivatives), and r, the radial coordinate at which the magnetic field needs to be evaluated. We recall that the magnetic field in the shell of the bubble given by Eq. (13) is defined only for r ∈ [R_min, R_max].

Our choice for the displacement field implies that a test particle initially at the origin (r⁰ = 0) is swept up to the position r = R_min and that the test particle initially at r⁰ = R_max remains at its location with r = R_max. With this model, the matter density and the magnetic field vanish inside the bubble (r < R_min) and remain unchanged outside it (r > R_max). Mass conservation from before the explosion to the present-time implies $$\rho (r)r^2 \text{d}r={\rho }^0\left(r^0\right)r^{0^2} \text{d}{r}^0,$$(19)

where ρ⁰ is the initial matter density. Assuming an initially homogeneous matter density and particularizing to our choice for Λ, the current matter density distribution is a function of r and reads as: $$\rho (r)=\{\begin{array}{ll}0\hfill & \text{ if }r<R_{\min }\hfill \\ {\rho }^0{\left(\frac{R_{\max }}{R_{\max }-R_{\min }}\right)}^3{\left(1-\frac{R_{\min }}{r}\right)}^2\hfill & \text{ if }r\in \left[R_{\min },R_{\max }\right]\hfill \\ {\rho }^0\hfill & \text{ if }r>R_{\max }.\hfill \end{array}$$(20)

Matter density varies with r inside the shell, and can also vary as a function of angular coordinates if R_min and R_max vary, that is, if the displacement vector field is not spherically symmetric.

For the sake of simplicity, we disregard any inhomogeneities or large-scale gradients in the initial matter density field. To deal with them, the initial density field and the displacement vector field would have to be modeled jointly, using the continuity equation, for example. We defer such a consideration to future work. Therefore, in its current implementation, although it enables us to take into account complicated shapes for bubble shells that might result from inhomogeneities (Kim & Ostriker 2015), our model best describes the regular component of the magnetic field in the shell of bubbles that result from a single explosion in a homogeneous medium, but possibly with anisotropic energy deposition.

Below, we look into the properties of the present-time magnetic field for the case of spherical bubbles with explosion centers at their origins (Sect. 2.3), and then for the case where spherical symmetry is lost (Sect. 2.4). In all cases, and as expected, the magnetic field in the shell of the bubble undergoes an amplification. We also demonstrate in Appendix A that our analytical model leads to magnetic flux conservation.

2.3 The case of spherical bubbles

In this subsection, we particularize our equations to the case where the displacement vector field is spherically symmetric. In this case, the angular derivatives of R_min and R_max vanish. This simplifies the expression for B_r in Eq. (13). It is then possible to write a simple expression for the amplification of the magnetic field in the shell of the bubble as due to the compression of the field lines. Substituting Eqs. (16) and (17) in Eq. (13) where we have ∇_tΛ = 0 for spherical bubble, we obtain $$\begin{array}{l}{\left(\frac{B}{B^0}\right)}^2={\left(\frac{r-R_{\min }}{r}\right)}^2{\left(\frac{R_{\max }}{R_{\max }-R_{\min }}\right)}^4\\ \times \left({\sin }^2\alpha +{\cos }^2\alpha {\left(\frac{r-R_{\min }}{r}\right)}^2\right),\end{array}$$(21)

where we introduced α, the angle between the direction of the initial magnetic field (B⁰) and the unit radial vector e_r pointing away from the explosion center. We have α = arccos(e_r · B⁰/|B⁰|).

In the triptych of Fig. 1, we show the crosscuts through an example of a spherical bubble where the amplification of the magnetic field can be visualized along with the projections of B in the respective planes. In this example, where B⁰ belongs to the XZ plane, the radial and zenith-angle dependence of the amplification can be appreciated. Equation (21) makes it clear that, for a fixed value of α, the strength of the magnetic field (B = |B|) increases from 0 at r = R_min to its maximum at r = R_max. We notice that this is a major difference as compared to models for the magnetic field in the shell of spherical bubbles constructed from geometrical arguments (e.g., van der Laan 1962; Wolleben 2007; Vidal et al. 2015; Korochkin et al. 2025), where the strength of the magnetic field is radially constant across the shell. Our model better matches results from magnetohydrodynamic (MHD) simulations of supernova explosions (e.g., Kim & Ostriker 2015; Maconi et al. 2023). In the outer surface of the bubble shell, the minimum value for B occurs at the magnetic poles, for α = 0 or 180°, and is such that B = B⁰. That is, there is no amplification of the field parallel to the initial magnetic field. This can also be observed in the middle panel of Fig. 1.

The dependence of the amplification of the magnetic field as a function of r and α are further illustrated in Fig. 2 for the same example spherical bubble with an inner radius of 100 pc and a shell thickness of 20 pc. The maximum amplification happens in the directions perpendicular to B⁰ , where the field lines are squeezed the most by the explosion. For α = 90° and r = R_max, the strength of the field at the current time is

Fig. 1 Crosscuts along the planes XY, XZ, and YZ through a spherical bubble with an inner radius of 100 pc and a thickness of 20 pc. The explosion center (black asterisk) is at the origin of the axes and at the center of the bubble. The direction of B0 is set by $\left({\theta }_{B^0},{\varphi }_{B^0}\right)=\left({60}^{°},0^{°}\right)$ in a spherical coordinate where θ = 0 points to positive Z and (θ, ϕ) = (0, 0) points to positive X. B0 belongs to the XZ plane. The background colors indicate the amplification of the total magnetic field in the plane (|B|/|B0|) and streamline the orientation of the projected B in the crosscut planes. The red arrow in the lower left corner of each panel indicates the projection of B0 in the planes. The red circle has a radius of B0 and indicates the maximum length of the arrow

Fig. 2 Dependence of the amplification factor of the magnetic field in the shell of a spherical bubble as a function of the radial distance (top) and as a function of the angle α between the radial unit vector and the direction of the initial magnetic field (bottom), for several values of α and r, respectively. In this example, Rmin = 100 pc and Rmax = 120 pc.

$$B=B^0\left(\frac{R_{\max }}{R_{\max }-R_{\min }}\right).$$(22)

Therefore, the maximum strength of the magnetic field depends on the size of the bubble and the thickness of its shell. The larger the bubble and the thinner its shell, the larger the amplification of the magnetic field strength.

By construction, the magnetic field inside the bubble (r < R_min) is set to zero while it is unchanged outside the bubble. Therefore, the amplification factor is respectively 0 and 1 in these regions. We note that, due to our choice for the displacement vector field, the final magnetic field is radially discontinuous at r = R_max. This is true for any bubble shape. This discontinuity could be resolved by choosing a continuous form for the radial dependence of Λ, which we will address in future work.

Finally, we note that our model predicts a relation between the amplitude of the magnetic field and the density within the bubble shell. The comparison of Eqs. (20) and (21) indicates that the B−ρ relation depends on the angle α. We have B ∝ ρ^1/2 at the magnetic equator, when α = 90°. The exponent slowly increases as α decreases and asymptotically converges to the value of 1 at the magnetic poles (α = 0° or 180°). As such, the explosion- induced compression of the matter density and amplification of the strength of the magnetic field imply a B–ρ relation that is close to B ∝ ρ^1/2 for most of the locations in the shell of the bubble. This relation corresponds to the case of lateral compression of elongated filaments of the ISM that are not parallel or perpendicular to their permeating magnetic field lines (e.g., Tritsis et al. 2015; Seta & Federrath 2022). It also agrees with results from MHD simulations of supernova explosions (e.g., Maconi et al. 2023).

2.4 The nonspherical cases

Our equations for the magnetic field in the thick shell of bubbles make it possible to consider cases where the spherical symmetry is lost. This may happen for spherical bubbles but with an off-centered explosion center or, more generally, to nonspheri- cal bubbles. The general case should occur in the ISM due to a certain degree of anisotropy in the initial velocity field related to the explosion of the progenitor supernova, to inhomogeneous initial density in which the supernova remnant develops (Kim & Ostriker 2015), or also to other feedback reactions from the compressed matter (e.g., Ferrière et al. 1991; Lopez et al. 2011; Orlando et al. 2016). To determine the magnetic field in such cases, it is necessary to estimate the angular derivatives of the inner and outer surfaces of the bubble shells. For this purpose, it is worth noticing that the orthoradial components of the gradient of the surfaces are also given by the orthoradial components of the normal vectors to the surfaces. Introducing n(θ, ϕ), the unit normal vector to a close surface R(θ, ϕ) of ℝ³ that surrounds the origin, we have: $${\nabla }_{t}R=e_r-\frac{n}{e_r\cdot n}$$(23)

in the direction e_r given by (θ, ϕ). Therefore, the two orthoradial components of the gradient of R in the direction e_r are related to the normal vector as follows: $${\nabla }_t\text{Λ}\cdot e_{\theta }=\frac{1}{r}\frac{\partial \text{Λ}}{\partial \theta }=-\frac{n\cdot e_{\theta }}{n\cdot e_r}$$(24) $${\nabla }_t\text{Λ}\cdot e_{\varphi }=\frac{1}{r\sin \theta }\frac{\partial \text{Λ}}{\partial \varphi }=-\frac{n\cdot e_{\varphi }}{n\cdot e_r}.$$(25)

If R_min and R_max can observationally be determined for a given bubble shell (e.g., see Appendix B), and if the explosion center is given, the magnetic field at a location (r, θ, ϕ) in the shell can be estimated through Eq. (13) using the normal vectors of R_min and R_max toward (θ, ϕ) and substituting the above equations in Eq. (18). It is interesting to note that, even if the angular derivative of R_min and R_max differ for a given direction, what happens as soon as the shell thickness varies, the orthoradial component of the magnetic field (as seen from the explosion center) will keep the same throughout the shell thickness and is fully determined by the orthoradial component of the initial magnetic field. This is because $B_{\varphi }/B_{\theta }=B_{\varphi }^0/B_{\theta }^0$, as derived from Eq. (13).

In Fig. 3 we show examples of magnetic fields obtained for cases where the spherical symmetry is lost. We consider the same spherical bubble as in Fig. 1, but with an off-centered explosion, and a more convoluted bubble shape for two possible positions of the explosion center. The comparison of these examples, also with the middle panel of Fig. 1, illustrates very well the leading dependence of the magnetic field amplification ratio with (i) the amount of swept-up matter (proportional to R_max as measured from the explosion center) and (ii) the thickness of the shell (inversely proportional to R_max – R_min). It is worth noting, however, that for nonspherical bubble shells, the gradient term $\left({\nabla }_t\text{Λ}\cdot B_t^0\right)$ in the expression for B_r alters the simple picture given by Eqs. (21) and (22). As a result, and as also clearly illustrated in Fig. 3, an off-centered explosion center may lead to significant asymmetry in the amplitude of the magnetic field. This asymmetry also depends on the geometry of the bubble shell.

Fig. 3 Crosscuts along the XZ plane through three bubble cases: a spherical bubble with an off-centered explosion center (left), a nonspherical bubble with the explosion center at the origin of the axes (middle), and that off-centered (right). The explosion center is marked with a black asterisk and the origin of the axes with a black dot. The bubble shape for the middle and right panels are the same. The direction of the initial magnetic field (also marked with the red arrow) is the same for every panel, and as in Fig. 1. B0 belongs to the XZ plane. Colors and streamlines follow the same convention as in Fig. 1.

3 Application to the Local Bubble shell and its contribution to polarization observables

Among all the bubbles and super-bubbles in the ISM, the Local Bubble is certainly a special case as the Solar system currently resides inside it. In this second part of the paper, we use our analytical model to carry out a phenomenological study with the magnetized shell of the Local Bubble as a study case. In particular, we focus on the impact of the shape of the shell and the choice of the effective explosion center on predictions of the Local Bubble shell’s contribution to sky maps of Faraday RM and synchrotron polarized emission. In this theoretical paper, we do not intend to model the full sky based on this model, nor incorporate it into models for the large-scale GMF and perform a full exploration of the parameter space. Such an endeavor is deferred to future work. Instead, we merely estimate the possible contribution of the Local Bubble’s thick shell to observables for specific study cases. To this end, and as explained below, we combine publicly available data on the shape of the Local Bubble from dust data, and reasonable values found in the literature for the free parameters to describe properties of the ISM and explosion. Therefore throughout this section we assume that the Local Bubble as it stands today can be described by a single, effective explosion that happened in a homogeneous medium permeated by a uniform magnetic field. This is a working assumption that only data can refute. We ignore possible effects from the turbulent component of the magnetic field and from inhomogeneities in the initial matter density. These are crude simplifications, the effects of which will need to be studied in future work.

To predict the contribution from the thick shell of the Local Bubble to the Faraday RM and polarized synchrotron emission using our analytical model (Eqs. (13) and (16)), the following ingredients are needed. Firstly, the strength and direction of the initial magnetic field (before the explosion) need to be set. This is done by fixing the three parameters B⁰ , $l_{B^0}$ and $b_{B^0}.l_{B^0}$ and $l_{B^0}$ are the Galactic longitude and latitude of the spherical coordinate system centered on the Sun. An initial magnetic field with $b_{B^0}={90}^{°}$ would point to the north Galactic pole and $\left(l_{B^0},b_{B^0}\right)=\left(0^{°},0^{°}\right)$ would point toward the Galactic center. The two angles give the initial directions of the magnetic field and B⁰ its amplitudes. Secondly, the shapes of the inner and outer surfaces of the shell need to be given as input. As explained in Appendix B, we choose to use shell surfaces derived by Pelgrims et al. (2020) from the 3D dust map of Lallement et al. (2019). Thirdly, the location of the effective explosion center needs to be set. This is done by setting its Cartesian coordinates (x_c, y_c, z_c) in the (Galactic) Heliocentric reference frame where the Galactic center (longitude 0) is toward positive x and the Galactic north pole is toward positive z. Using our analytical model, the three above inputs (six parameters plus inner and outer surfaces) make it possible to estimate the magnetic field anywhere in the presenttime shell of the Local Bubble. Finally, to produce maps of the Faraday RM and of the synchrotron Stokes parameters Q and U, we need to adopt models for the thermal electron density, for the cosmic-ray electron (CRE) density, and for the energy-spectrum of the CRE. Our choices are described next.

3.1 Study cases

As evidenced in Sect. 2, the present-time magnetic field in the shell of the Local Bubble may severely depend on the specific geometry of the shell (shape, size, and thickness) and the position of the explosion center. We investigate the possible differences that specific choices imply on polarization observables by constructing six scenarios.

The first type consists of spherical bubble shells (with constant thickness). Their names start with SC. For the second type, the bubble shells have their inner surface directly corresponding to the model derived by Pelgrims et al. (2020) from 3D dust data and obtained with a maximum multipole l_max = 6 (i.e., $r_{\text{inner }}^{{\ell }_{\max }=6}$). Their names start with a D. For this second type, we consider the case where the shell thickness (as measured radially from a certain position) is constant (named DC) and the case where it varies (named DD). For this latter case, we adopt the outer surface of the shell $\left(r_{\text{outer }}^{{\ell }_{\max }=6}\right)$ as determined in Appendix B, also from the same 3D dust data.

To fix the center and the inner radius of the sphere for the spherical bubble-shell scenarios, we choose to fit a sphere to the 3D data points drawn from $r_{\text{inner }}^{{\ell }_{\max }=6}$. As such, our spherical model for the Local Bubble is also derived from 3D dust data. In the Heliocentric Cartesian coordinate system, the center of the sphere is found at (x_sph, y_sph, z_sph) = (–24.8, –32.6, –23.3) pc and its radius is R_sph = 216.7 pc. For the constant thickness of the shell, we adopt the value of ∆ = 35 pc, in agreement with the model of Yao et al. (2017), and that is not too far from the best-fit value of 30 pc obtained in Korochkin et al. (2025). For the case named DC we choose (x_sph, y_sph, z_sph) as the point from which the shell thickness, as measured radially, is constant.

Pelgrims et al. (2020) used $r_{\text{inner }}^{{\ell }_{\max }=6}$ to model the high-Galactic latitudes (|b| ≥ 60°) of the dust polarized emission measured by the Planck satellite at 353 GHz. They obtained constraints on the position of the center of the effective explosion that has led to the Local Bubble, and on the orientation of the initial magnetic field. Because dust-polarized emission is not directly sensitive to the amplitude of the magnetic field, their solution for the position of the explosion center is degenerate along a line parallel to the initial direction of the magnetic field. We adopt the solution that they found maximizes the likelihood (obtained using $r_{\text{inner }}^{{\ell }_{\max }=6}$), and that approximately reads as $\left(x_c^{\text{P}20},y_c^{\text{P}20},{ɀ}_c^{\text{P}20}\right)=(23,-34,-122)$ pc and $\left(l_0^{\text{P}20},b_0^{\text{P}20}\right)=\left({73}^{°},{17}^{°}\right)$.

We notice that the explosion center found in Pelgrims et al. (2020) is significantly different than the center of the sphere that best fits $r_{\text{inner }}^{{\ell }_{\max }=6}$. Therefore, we consider both locations as possible explosion centers to build our study cases. The third letter of our scenarios’ names indicates this choice. The letter O corresponds to explosion center at (x_sph, y_sph, z_sph) and letter A corresponds to explosion center at $\left(x_c^{\text{P}20},y_c^{\text{P}20},{ɀ}_c^{\text{P}20}\right)$.

In total, we are thus considering six scenarios to compute the present-time magnetic field in the thick shell of the Local Bubble. The definition of these scenarios, with our choices for their parameter values, are summarized in Table 1.

To fully determine the magnetic field in the shell of the Local Bubble for our scenarios, we need to set the value for the strength of the initial magnetic field, before the explosion happened. In the ensemble of models for the large-scale Galactic magnetic field presented by Unger & Farrar (2024), the coherent component has a field strength that ranges from 0.2 to 0.5 µG at the Sun location. However, due to the turbulent component of the magnetic field, the actual strength for the initial magnetic field may have been much larger. We adopt a value of 3 µG for the strength of the initial field (B⁰). This is on the lower limit of the turbulent field strength of 3 to 6 µG found for the disk by Planck Collaboration Int. XLII (2016) who updated the model of Jansson & Farrar (2012b). On the other hand, Korochkin et al. (2025) found a value of 3.5 µG for the strength of the GMF in the Solar neighborhood, which is close to our choice for B⁰ . We emphasize on the fact that, in our simple model, the strength of the initial magnetic field acts only as a scaling factor. A scan over possible values for B⁰ in the range from zero to 5 µG is also performed in Sect. 3.4.

3.2 Faraday RM and synchrotron Stokes parameters

Free (thermal) electrons of the magnetized ISM induce a net rotation of the polarization plane of an electromagnetic wave passing through the medium before it reaches the observer. This Faraday rotation is proportional to the square of the wavelength and is characterized by the RM. The RM results from the line- of-sight integration, from the observer to the source at distance d, of the product of the thermal electron and the line-of-sight component of the magnetic field. The RM is computed as $$\text{RM}\propto -{\displaystyle {\int }_0^d \text{d}}r{n}_e(r)B_{\Vert }(r),$$(26)

where n_e is the number density of thermal electrons and B_|| = B ⋅ e_r is the line-of-sight component of the magnetic field, where er is the radial unit vector pointing outward from the origin of the spherical coordinate system centered on the observer. To express the RM in rad/m² , and if the magnetic field is expressed in microgauss, the distance in parsec, and the thermal electron in cm^–3, the proportionality factor takes the value of ≈0.81 (rad/m²)(cm³/pc)(µG)^–1.

To predict the contribution of the Local Bubble shell to the RM sky map, we thus need to adopt a model for the thermal electron density. Thermal electron density models of Cordes & Lazio (2002) and Yao et al. (2017) predict a local value for n_e of 0.011 cm^–3 and 0.016 cm^–3, respectively. These values were obtained for the diffuse ISM, without the density depletion within the Local Bubble as implemented in these models. The mid-plane density of the electron density in the local Galactic disk as determined by Ocker et al. (2020) is n_e = 0.015 cm^–3. It is within the range of the two above models and we adopt it as our default value for the initial, homogeneous thermal electron density $\left(n_e^0\right)$, that is, for the thermal electron density before the effective supernova explosion that has led to the Local Bubble. According to our simple model for the formation of the Local Bubble, the thermal electron density distribution is modified by the explosion. We assume that the thermal electron density follows Eq. (20) and that the fraction of ionized to total gas remains constant so that the shell properties as inferred from 3D dust map apply for the thermal electron density. We notice that in this framework, changing the explosion center inside a given geometry for a bubble shell leads to different thermal electron densities in the shell.

Relativistic electrons in the ISM spiral about the ambient magnetic field lines and lose energy through synchrotron emission that is linearly polarized perpendicular to the orthoradial component of the magnetic field (Eq. (15)). Relativistic electrons with energies on the order of 10 GeV are responsible for the Galactic synchrotron emission observed in the frequency range of 10 to 30 GHz, well covered by cosmic microwave background experiments, and frequencies at which Faraday rotation is negligible.

We consider the simplified assumption that CRE follows approximately a power-law energy distribution of the form: $$N(E)\text{\hspace{0.17em}d}E=\kappa E^{-p} \text{d}E=n_{10}{\left(\frac{E}{10\text{\hspace{0.17em}GeV}}\right)}^{-p} \text{d}E,$$(27)

where p is the electron spectral index, and where we adopt a normalization at 10 GeV energy. In this case, the energy spectrum of CRE is related to the cosmic-ray flux by: $$n_{10}{\left(\frac{E}{10\text{\hspace{0.17em}GeV}}\right)}^{-p}=\frac{4\pi }{c}{\text{Φ}}_{10}{\left(\frac{E}{10\text{\hspace{0.17em}GeV}}\right)}^{-p},$$(28)

where Φ₁₀ = 10^–0.55 Gev^–1 m^–2 s^–1 sr^–1 at 10 GeV. Therefore, we have n₁₀ = 1.2 × 10^–23cm^–3eV^–1 and κ ≈ 3.03 × 10^–25 J² m^–3 for p = 3. These numbers were obtained by Unger & Farrar (2024) from the flux in the local ISM (i.e., outside of the heliosphere) as derived from AMS2 data (Aguilar et al. 2014, 2015, 2019). Following Longair (2011), assuming p = 3 and a uniform CRE density, the simulated Stokes parameters Q and U for the linearly polarized synchrotron emission at a given frequency ν can be written in terms of line-of-sight integration over solid angle of the plane-of-sky components of the magnetic field as $$\begin{array}{l}{Q}_v=-\frac{S_v}{4\pi }{\displaystyle {\int }_0^{\infty }\text{d}}r\left(\left(B_{\theta }{(r)}^2-B_{\varphi }{(r)}^2\right)\right)\\ U_v=-\frac{S_v}{4\pi }{\displaystyle {\int }_0^{\infty }\text{d}}r2B_{\theta }(r)B_{\varphi }(r),\end{array}$$(29)

and can be converted to antenna temperature for comparison with data. For observation at 30 GHz, and if we express the line-of-sight integration in parsec and the magnetic field in microgauss, the amplitude S_ν in Eq. (29) is ≈2.8 × 10^–27 J/(pc m² µG²). The “–” signs in the above expressions are such that the Stokes parameters are given in the HEALPix convention (positive Q to the south, negative Q to the east, positive U to the southeast, negative U to the northeast), also referred to as COSMO convention (Górski et al. 2005).

Although cosmic rays are not passive during the formation of a supernova remanent (e.g., van der Laan 1962), and may even trigger deformation of the bubble shell due to pressure balance (e.g., Ferrière & Zweibel 1991), we do not consider these effects in this exploratory study. Instead, we assume that, at the present time, the CRE are distributed uniformly as described above. Further work will be carried out to go beyond this simplifying assumption. If CRE density has been enhanced in the shell of the Local Bubble, our calculations provide a lower limit to the possible contribution to the synchrotron polarized emission.

3.3 Results

The full-sky maps of the predicted contribution of our bubbleshell scenarios to the Faraday RM and Q and U Stokes parameters of the polarized synchrotron emission are shown in Fig. 4. The maps were produced using an HEALPix angular tesselation with N_side = 64 (Górski et al. 2005). For the radial sampling, we use 100 grid points through the thick shell of the Local Bubble as seen from the Sun. The line-of-sight integrations (Eqs. (26) and (29)) were performed using the simple midpoint rule and using the parameter values given in the previous subsections and summarized in Tables 1 and 2. We checked that increasing the radial sampling does not affect our result.

As seen in Fig. 4, and as expected, the convoluted shape of the Local Bubble shell leads to a certain number of anisotropies in the maps, the details of which depend on both, the choice of the explosion center and the thickness of the shell (in addition to its shape). The shape and exact position of these anisotropies in the sky also depend on the orientation of the initial magnetic field (not shown in the figure). The substantial north-south asymmetry seen in the maps of the DCA and DDA cases as compared to DCO and DDO cases is due to the fact that in A cases, the explosion center is located significantly more below the Galactic plane (–122 pc as compared to –23.3 pc). As a consequence, the amount of matter and magnetic field lines that have been swept up to the northern part of the Bubble is much larger than to the south. This is reminiscent of what was discussed in Sect. 2 and could, to some extent, be at the origin of the asymmetry observed in RM data (Dickey et al. 2022).

To compare the model predictions between them and also compare them with observations, we show longitude profiles for several stripes of constant latitudes in Fig. 5. We consider latitude stripes of the width of 10° with centers at 60°, 70°, and 80°, both in the northern and southern hemispheres. The longitude profiles are constructed by taking the average of the model predictions, or observations, in longitude bins of 10° width. The errorbars corresponding to the data show the standard deviation of the observation in each bin. It is indicative of the scatter in the data that primarily results from small-scale fluctuations in the magnetized ISM (e.g., Jaffe et al. 2010; Jansson & Farrar 2012a; Pelgrims et al. 2021; Unger & Farrar 2024; Korochkin et al. 2025), plus the source contribution in case of Faraday RM. For this comparison, we use the RM data from Unger & Farrar (2024) at N_side = 32, and the synchrotron Q and U data at 30 GHz from Planck processed data by the Commander component separation method (Planck Collaboration IV 2020), smoothed at a resolution of 1.2°, and downgraded at N_side = 64. We do not mask the data to build these longitude profiles. For example the North Polar Spur is well seen in the region with l ≲ 90° and 55° < b < 75° in both Q and U data.

In Fig. 5 it is seen that the different scenarios lead to variations of the RM signal whose number and amplitude are model-dependent. They depend on the shape of the shell and on its thickness. Overall, the contribution from the Local Bubble shell to the RM sky is smaller (in amplitude) than the observed signal. This indicates that the Local Bubble shell is on average only a subdominant component of the RM integrated over the full path length through our Galaxy. This result is at odds with that of Reissl et al. (2023) who found, based on MHD simulations of Milky Way-like galaxies, that ISM bubble shells surrounding an observer contribute significantly to the RM signal at high latitudes. This may indicates that we are either underestimating the value of $\left(n_e^0\right)$ or missing a contribution to the line-of-sight component of the magnetic field in the shell that could come from the compressed turbulent component. However, it is interesting to note that in some parts of the sky, such as at b > 75° for l < 90°, the signal from the shell of the Local Bubble that is predicted for some of our scenarios reaches the amplitude in the data and, furthermore, that the RM data are not very constraining at these latitudes. Overall, the significance of the discrepancy between our model predictions and the data is relatively low, as further discussed in Sect. 3.4.

The shape and thickness of the Local Bubble shell, as well as the choice for the location of the explosion center, also induce changes in the predictions for the Q and U Stokes parameters of the synchrotron polarized emission. Overall, the different scenarios lead to the same phase for the sinusoidal variation of the signal with longitude. This is because we assume the same direction for the initial magnetic field. The different scenarios, however, predict variations in the amplitude of this sinusoidal signal and even departure from a perfect sine function. These differences primarily depend on the assumed shape and thickness of the shell.

Although the contribution from the Local Bubble shell is only subdominant for the RM at those latitudes, this is not the case for the synchrotron Q and U, as seen in Fig. 5. Surprisingly, the contribution from the Local Bubble shell to the synchrotron sky is very significant for |b| ≳ 55°. Of course, alone the shell of the Local Bubble cannot account for the entirety of the signal and its variations. However, it is very striking that our simple model with our a priori choice for the model parameters leads to a shell contribution to the Q and U sky that has the same overall amplitude and phase as that of the observation. It should be emphasized that the models do not result from a fit to RM or synchrotron data. The result of a fit would likely improve the agreement seen in Fig. 5, especially for Q and U. Further analysis will be needed to better understand why the agreements between data and model predictions are different for Faraday RM than for synchrotron polarized emission. Such an analysis might require the inclusion of the effects of the turbulent component of the magnetic field and inhomogeneities, and will require the full exploration of the parameter space. This is beyond the scope of this exploratory study.

Fig. 4 Mollweide projection of the full-sky maps of the contributions from the shell of the Local Bubble to the RM, Q, and U signal (from left to right) as predicted for our different scenarios (rows). Maps are given in Galactic coordinates with the Galactic center at the center of the maps, Galactic longitude increases to the left and the Galactic north pole is at the top. Color scales are shared in columns. They range from -10 to 10 rad/m2 for the RM maps and from –10 to 10 µK for Q and U maps.

Fig. 5 Longitude profiles for constant latitude stripes of 10° width for |b| ∈ [55°, 85°], as described in the text. Predictions for the different scenarios for the Local Bubble magnetic field are shown as continuous lines according to the legend, for the RM, Q, and U signal (from left to right). For reference, observational data is also shown as gray points with symmetric errorbars.

3.4 Discussion

Although we do not attempt a full adjustment of the model to the data, we proceed to a better quantification of the comparison between the maps from the different scenarios of the Local Bubble shell and the data. We consider the same data as used by Unger & Farrar (2024) to constrain models of the large-scale regular component of the Galactic magnetic field. We used their RM, Q, and U maps at the N_side = 16 resolution, along with their uncertainties, and their masks. For our comparison, we choose to keep only all unmasked pixels with Galactic latitude |b| > 45° to avoid sky regions where the signal is clearly dominated by other Galactic components than the wall of the Local Bubble, such as the Galactic disk. We are left with 751 pixels for the Q and U maps and with 819 pixels for the RM map. The masked data are shown at the top row of Fig. 6. Accordingly, we downgrade the modeled maps shown in Fig. 4 from N_side = 64 to N_side = 16 and apply the same mask than to the data. If $d_X^i,m_X^i,$ and ${\sigma }_X^i$ are the data, model and uncertainty for the observable X in pixel i, where X is either RM, Q, or U, we then compute the residual $\left(d_X^i-m_X^i\right)$ and pull $\left({\chi }_X^i:=\left(d_X^i-m_X^i\right)/{\sigma }_X^i\right)$ for each pixel. Using these values, we can compute the reduced χ² of each map and for each scenario. The values are given in Table 3, along with a total value that results from the combination of RM, Q, and U maps. The maps of the model and pull for the scenario SCA are shown in Fig. 6 at the second and third rows, respectively.

The qualitative match found between the data and the six scenarios in Fig. 5 is confirmed. Among them, the SCA scenario provides the best performances and the DCA scenario the worst, with a total reduced χ² of 0.95 and 1.7, respectively. These χ²/ndf are obtained without any minimization. In addition, the low values found for the reduced χ² are due to the relatively large uncertainties assigned to the data at those latitudes, in particular for the Faraday RM as discussed below. We notice that the pull maps of Q and U shown in Fig. 6 show very well-structured patterns that are unlikely to be the result of random noise. This is evidence that other large-scale components need to be included to model the data at those latitudes, such as a halo field for example.

To investigate further the relevance of the contribution of the shell of the Local Bubble to the RM, Q, and U maps, we vary the strength of the initial magnetic field (B⁰) from 0 µG, where the Local Bubble shell has no contribution to the observable, to a value of 5 µG. For each value of B⁰ , we compute the reduced χ² for RM separately and Q and U combined, and also consider the northern and southern hemispheres at a time. If the data has sufficient constraining power, and if a model contributes substantially to the data with a relatively good match, we expect to observe a well-defined minimum in the curve of reduced χ² versus B⁰ . If, instead, the model contribution is negligible for all B⁰ values, we expect to see a relatively flat curve.

The results for the six scenarios under consideration are shown in Fig. 7. On the one hand, it is seen that Faraday rotation contributes only marginally to the total reduced χ². For all tested values of B⁰, the reduced ${\chi }_{\text{RM}}^2$ varies only marginally and remains close to the value where the contribution from the Local Bubble shell is assumed to be zero (i.e., when B⁰ = 0). This value is found to be ≈0.9, showing that the S/N of RM data for |b| > 45° is low on average. This demonstrates that the constraining power of the RM data at those latitudes is rather weak. On the other hand, for all six scenarios, we observe a clear minimum in the curves corresponding to Q and U, and the reached minima are significantly lower than in the case of no contribution from the shell of the Local Bubble. For all scenarios, the minima are found for a value of B⁰ in the range of 2 to 4 µG. For the cases with the explosion center at the position O (see Table 1), the curves for the northern and southern hemispheres are very similar. However, for the cases with the explosion center at the position A, the agreement between the model and the data is better in the northern hemisphere (dashed line) than in the southern hemisphere (dotted line). Indeed, the minima are more pronounced in the northern part than in the southern part of the sky. There is also a tension between the favored value for B⁰ from both hemispheres. To understand better these features, the magnetic field in the shell of the Local Bubble needs to be included in global models for the regular component of the GMF, and a full exploration of the parameter space is needed. This task goes well beyond the scope of this paper. Overall, this latter test confirms that the shell of the Local Bubble likely makes a significant contribution to the observed Q and U Stokes parameters of the Galactic synchrotron emission whereas it only makes a marginal contribution to the Faraday RM integrated over the full path length across the Galaxy. This is obtained for Galactic latitudes |b| > 45°, considering the masks of Unger & Farrar (2024) that, in particular, masks out bright features of the synchrotron sky, and without having explored the parameter space of our model other than the value of B⁰ . Our latest test also enabled us to test explicitly the robustness of our conclusion against the specific choice for the value of the strength of the initial magnetic field. In fact, for a fixed bubble-shell geometry, fixed explosion center, and fixed orientation of the initial magnetic field, and according to our model assumptions, the contribution of the Local Bubble shell to the RM scales with $n_e^0{B}^0$, while its contribution to the Q and U scales with n₁₀ (B⁰)². Consequently, one way to increase the RM contribution of the Local Bubble shell, and thus better match the RM data without deteriorating the qualitative agreement observed for the synchrotron polarized emission, would be to increase substantially the ratio $n_e^0/n_{10}$ as compared to the one we used in this work. Either $n_e^0$ needs to be higher, or n₁₀ to be smaller. If n₁₀ overestimates the true density of CRE in the shell of the Local Bubble, then a larger value for the magnetic field would be needed. This would allow for an increased contribution to the RM.

To conclude this discussion section, we emphasize that, although it improves on the current state-of-the-art, our model remains simple and is based on several assumptions that may not apply satisfactorily to the case of the Local Bubble, and whose effects are not all known. For instance, variations in both the magnetic field and the initial matter density are completely disregarded in our model. It could be argued that fluctuations that are then stretched and squeezed during the bubble formation could change the balance between the predicted amplitudes of the Faraday RM and synchrotron polarized emission (e.g., Beck et al. 2003). This could possibly reconcile our results with those from MHD simulations of Milky Way-like galaxies where the observer is located inside super-bubbles (Reissl et al. 2023), or not (Pakmor et al. 2018). Addressing such open questions is beyond the scope of this paper and will be the subject of follow-up studies.

Fig. 6 Orthographic projection of the same data used in Unger & Farrar (2024) (top), the model predictions for the SCA scenario (middle), and pulls (bottom). The gray area is masked out as explained in the text. The north pole is at the center of the left disk and the south pole to the right. The longitude zero is marked by the vertical thick lines starting from the poles upward. Longitude increases clockwise on the left and counterclockwise on the right.

Fig. 7 Contribution to the reduced χ2 of the different observables as a function of the strength of the initial magnetic field for the six scenarios. The contribution from Q and U are combined. The contributions from the northern and southern hemispheres are also shown with dashed and dotted lines, respectively.

4 Summary and concluding remarks

An analytical prescription for the coherent magnetic field in the thick shell of ISM (super-)bubbles has been derived. We constructed our analytical model starting from the induction equation without magnetic diffusion and the basic assumption that ISM bubbles result from a single explosion that induces a purely radial displacement of the surrounding gas, but not necessarily with spherical symmetry. Our model assumes uniform initial conditions. Relying on the frozen-in approximation, we obtained equations for the magnetic field components expressed in spherical coordinates centered on the explosion center. The solution makes it possible to estimate the magnetic vector field anywhere in the thick shell of the bubble. By construction, the solution leads to a magnetic field that is divergence-free. To study and illustrate the properties of our analytical solution, we further proposed a linear model for the radial displacement vector field that relates the initial locations of matter particles to their present locations in the thick shell of the bubble. For the bubble cases with spherical symmetry, we demonstrated that our analytical model leads to the general expectations in terms of geometry for the magnetic field lines and amplification. Our solution also makes it possible to study other cases where the spherical symmetry is broken. We showed that strong anisotropies in magnetic field amplification can occur, not only due to the magnitude of the radial displacement but also due to the circumvolutions in the shell geometry.

Assuming that the inner and outer radius of ISM bubbles can be extracted from observation, our analytical prescription makes it possible to model the magnetic field in their shell, with some limitations and caveats. For example, due to the choice of the radially discontinuous, linear model for the radial displacement field, our magnetic field undergoes a sharp discontinuity at the outer radius of the bubble shell. This discontinuity could be overcome by adopting a continuous model for the radial displacement and will be addressed in future works. The main limitation of our model, in its general form or with the linear displacement field, is that it does not explicitly account for the dynamics of the bubble formation, but solely depends on the initial and final (presenttime) conditions. While this is the main limitation from the point of view of the physical modeling of ISM bubbles, we emphasize that it is also the main asset of our model for all practical uses.

We further used our model to explore the possible contribution of the shell of the Local Bubble to the extragalactic Faraday RM and synchrotron polarization data at 30 GHz, which our model makes it possible to do since it predicts the amplitude of the magnetic field in the thick shell. For this study, we relied on 3D geometrical models for the shell of the Local Bubble informed from a 3D dust density map and adopted published values for all free parameters in our model. We considered six scenarios with varying complexity for the bubble shell, considered either constant or varying shell thickness, and explored two possible explosion centers. We found that all these geometrical parameters lead to (sometimes large) variations in the maps of observables. For the six scenarios, we found that the shell of the Local Bubble may be a substantial contributor to the synchrotron data at high Galactic latitudes (|b| > 45°), whereas its contribution to the integrated Faraday RM is generally small. In that respect, our study reinforces the recent result of Korochkin et al. (2025) that indicates that, at high Galactic latitudes, the shell of the Local Bubble is an important foreground to the magnetized Galaxy and that its contribution to the synchrotron polarization data is significant while subdominant for the RM data. This comes with an improved model for the magnetic field in the shell of the Local Bubble and an improved model for the shell itself. In addition, we showed that the specific geometry of the shell and the location of the explosion center inside the bubble are crucial parameters that may lead to very different predictions. A spherical model can be used to approximate the overall patterns in the RM, Q, and U maps, but only the solution for the irregular Bubble can lead to detailed fine-structures seen in observation. When applied to the case of the Local Bubble, the main limitations of our model possibly come from the fact that we assume the Local Bubble to result from a single explosion that happened in a homogeneous medium permeated by a uniform magnetic field. Given the size of the Local Bubble, these are crude assumptions that will have to be move past in the future.

Our studies underline the need to include the Local Bubble and its magnetized shell in models of the regular component of the GMF. In fact, focusing on the high Galactic latitudes, we find that the predictions from the Local Bubble shell have signal amplitude comparable to the predictions from the halo components of the regular Galactic magnetic field used in Unger & Farrar (2024) and shown in their Fig. 10. Only a combined fit to the data, including the shell of the Local Bubble into the description of the large-scale GMF, with a full exploration of the parameter space will make it possible to quantitatively judge on their respective contribution. We will address such an endeavor in future work and explore the possible impact on the studies of the propagation of ultra-high energy cosmic rays across the Galaxy. In this context, the analytical expression for the magnetic field in the Local Bubble shell obtained in this paper will prove very useful, as it enables us to determine the magnetic field with only a handful number of parameters. Of course, this model relies on several working assumptions and does not explicitly account for all the physics involved in a supernova explosion such as pressure balance, feedback from ISM matter onto the development of the shell (e.g., van der Laan 1962; Ferrière et al. 1991; Ferrière & Zweibel 1991), and other phase transitions that may take place over the time in such peculiar environments. However, our model is general enough to effectively accommodate a diverse range of complexities and to be useful for better understanding magnetized bubbles and super-bubbles of the ISM, also other than the Local Bubble itself. Modeling magnetized supernova remnants may also prove useful in deciding between competing large-scale GMF models (West et al. 2016), a context in which our analytical model could present advantages.

Acknowledgements

We are grateful to the anonymous referee who provided us with a constructive report that helped us improve our manuscript. The authors would like to thank Katia Ferrière, Alexander Korochkin, and Peter Tinyakov for helpful discussions. VP acknowledges funding from a Marie Curie Action of the European Union (grant agreement No. 101107047).

Appendix A Magnetic flux conservation

In this appendix, we demonstrate that our analytical model for the magnetic field in thick shells of ISM bubbles with the linear displacement field leads to the conservation of the magnetic flux, before and after the explosion. We use the same notation and convention as in Sect. 2.

We start by using the case of a spherical bubble. The magnetic flux through a plane with unit normal vector u is given by Φ = ∫ dA B ⋅ u. To demonstrate that Φ = Φ⁰, we consider a plane that contains the explosion center and with a normal vector that makes an angle γ with B⁰ . To demonstrate that the flux is conserved we need just to integrate B ⋅ u over the surface delimited by the outer radius of the spherical bubble. The field is unchanged outside it.

Without loss of generality, we consider a spherical coordinate system centered on the explosion center and with the north pole pointing parallel to u. The measure of the flux thus happens in the equatorial plane of this coordinate system, where u ⋅ e_θ = −1. By construction, we have B⁰_θ = B⁰ ⋅ e_θ = −B⁰ cos γ. The initial flux through the surface is therefore given by Φ0 = πR_max² B⁰ cos γ.

To estimate the flux of the present-time magnetic field through the same surface, we need to integrate B ⋅ u. We have $$\text{Φ}={\displaystyle \int \text{d}}AB\cdot u$$(A.1) $$=-{\displaystyle \int \text{d}}A{B}_{\theta },$$(A.2)

where, for a spherical bubble and with the linear displacement model, B_θ (r) is given by $$B_{\theta }={\left(\frac{R_{\max }}{R_{\max }-R_{\min }}\right)}^2\left(\frac{r-R_{\min }}{r}\right)B_{\theta }^0$$(A.3)

in the thick shell, and is zero for r < R_min . Performing the integral, we have $$\text{Φ}=-{\displaystyle \int \text{d}}A{B}_{\theta }$$(A.4) $$=-{\displaystyle {\int }_0^{2\pi } \text{d}}\varphi {\displaystyle {\int }_{R_{\min }}^{R_{\max }}\text{d}}rr{B}_{\theta }(r)$$(A.5) $$=-2\pi B_{\theta }^0{\left(\frac{R_{\max }}{R_{\max }-R_{\min }}\right)}^2{\displaystyle {\int }_{R_{\min }}^{R_{\max }}\left(r-R_{\min }\right)}\text{\hspace{0.17em}d}r$$(A.6) $$=-\pi R_{{\max }^2}{B}_{\theta }^0$$(A.7) $$={\text{Φ}}^0.$$(A.8)

This shows that the flux of the magnetic field before and after the explosion is conserved in any plane containing the explosion center. This demonstration for the magnetic flux conservation can be generalized to the case of a nonspherical bubble shell. In this case, one needs to take into account the azimuth-angle dependence of the outer and inner radius of the shell in the equatorial plane. In the general case, the flux is given by the following: $$\text{Φ}=\frac{B^0\cos \gamma }{2}{\displaystyle {\int }_0^{2\pi } \text{d}}\varphi R_{\max }(\varphi ),$$(A.9)

where γ is the angle between the initial magnetic field direction and the normal to the plane in which magnetic flux is measured. This is exactly the expression for Φ⁰ for such a case.

Appendix B The thick shell of the Local Bubble

In the last decade, precise measurements of stellar parallaxes from the Gaia space mission (Gaia Collaboration 2016) and stellar reddening has enabled to map the dust density distribution in three dimensions (e.g., Green et al. 2019; Lallement et al. 2019; Leike & Enßlin 2019; Lallement et al. 2022; Vergely et al. 2022; Edenhofer et al. 2024). The advent of such 3D maps enables extracting the shell of the Local Bubble as a region of higher dust density and characterizing its geometry directly from observation. Pelgrims et al. (2020) were the first to extract the shape of the inner surface of the Local Bubble shell based on 3D dust maps. They considered the maps of Lallement et al. (2019) and Leike & Enßlin (2019), which were the most recent maps at the time and enabled a full-sky study to be carried out.

Drawing sightlines in the data cube outward from the Sun, they constructed radial profiles of the differential of the dust extinction δA_v ≡ dA_v /dr (with r the distance to the Sun). For each sightline separately, they defined the inner radius of the Local Bubble shell as the distance to the first inflection point where the profile transitions from convex to concave. Having obtained this measure for all sightlines of an HEALPix angular tesselation with N_side = 128, they then derived smoothed versions of the inner surface of the shell using spherical harmonic decompositions of the surface with limiting maximum multipole (ℓ_max). Pelgrims et al. (2020) used these smoothed models of the inner surface of the Local Bubble shell to successfully constrain a model for the smooth component of the magnetic field in the shell, which was assumed to be very thin, through a fit of the dust polarized emission at high-Galactic latitudes. They found that the inner surface extracted from the 3D map of Lallement et al. (2019) leads to more convincing results.

Although there have been new 3D maps of the dust density and also more recent characterization of the shape of the shell of the Local Bubble based on these maps (e.g., O’Neill et al. 2024b), we choose to use the model of Pelgrims et al. (2020) obtained with the maximum multipole ℓ_max = 6 for the inner surface of the Local Bubble¹. This choice is motivated by the smooth, closed nature of the surface, which makes it particularly suited to our purposes. We refer to this surface as $r_{\text{inner }}^{{\ell }_{\max }=6}$, both in this Appendix and in the main text. $r_{\text{inner }}^{{\ell }_{\max }=6}$ consists of a sky map of distances from the Sun to the inner edge of the shell of the Local Bubble.

To describe the outer surface of the Local Bubble shell, we consider two cases. For the first case, we assume that the thickness of the shell is constant as seen from a given origin located inside the bubble (see Sect. 3.1). Therefore, the outer surface of the shell is fully determined by $r_{\text{inner }}^{{\ell }_{\max }=6}$, the adopted origin, and a constant thickness ∆ that we fix to 35 pc, in agreement with the model of Yao et al. (2017) for the thermal electron density that includes a model for the Local Bubble. A value of 35 pc for the shell thickness gives a ratio of thickness to radius on the order of 10 to 30% that agrees with what is observed and predicted for ISM super-bubbles (e.g., Krause & Diehl 2014). However, and to mention only the most obvious, given the stratified and inhomogeneous environment in which the Local Bubble has formed and expanded, it is very unlikely that the thickness of its shell is constant. So, to consider the case of variable thickness, we use an approach very similar to that used to calculate $r_{\text{inner }}^{{\ell }_{\max }=6}$ in (Pelgrims et al. 2020) and derive an estimate of the outer surface of the Local Bubble shell from 3D dust data. We use the 3D dust map of Lallement et al. (2019) and study the radial profiles of the differential dust extinction. Instead of considering the first inflection point where the profile goes from convex to concave, we consider the second inflection point where, after reaching a local maximum, the profile transitions from concave to convex. We repeat the process for all sightlines of an HEALPix map with N_side = 128. Because of the limited extent of the 3D map, the second inflection point is either not found before the edge of the map or found very close to it for about 20% of the sightlines (Pelgrims et al. 2020). Therefore, to model the outer surface of the shell, we mask those problematic sightlines and fit the unmasked data by adjusting a surface that is fully described by the spherical harmonic coefficients up to a maximum multipole ℓ_max. For consistency with the characterization of the inner surface, we use ℓ_max = 6 and define the outer surface of the shell of the Local Bubble, which we denote $r_{\text{outer }}^{{\ell }_{\max }=6}$. Similar to $r_{\text{inner }}^{{\ell }_{\max }=6},r_{\text{outer }}^{{\ell }_{\max }=6}$ con- sists of a sky map of distances from the Sun to the outer edge of the shell of the Local Bubble. The maps of $r_{\text{inner }}^{{\ell }_{\max }=6}$ and $r_{\text{outer }}^{{\ell }_{\max }=6}$ are shown in the left and right panels of Fig. B.1, respectively. In Fig. B.2, we show three crosscuts of the Solar neighborhood according to the 3D dust extinction map of Lallement et al. (2019). This map gives the differential extinction (δA_v). These crosscuts show the XY, XZ, and YZ planes, of the Heliocentric Cartesian coordinate system where the Galactic center is located at X ≈ 8 kpc. In this triptych, the blue and orange points show the locations of the first and second inflection points measured on the radial profiles δA_v(r) in the respective plane, when they are well-defined. The continuous lines correspond to the intersection of our models for $r_{\text{inner }}^{{\ell }_{\max }=6}\text{and\hspace{0.17em}}{r}_{\text{outer }}^{{\ell }_{\max }=6}$ with the different planes, respectively.

Fig. B.1 Sky maps of the modeled distances to the inner $\left(r_{\text{inner}}^{{\ell }_{\max }=6},\text{left}\right)$ and outer $\left(r_{\text{outer}}^{{\ell }_{\text{max}}=6},\text{right}\right)$ surfaces of the shell of the Local Bubble as measured from the Sun. Projection and conventions are the same as in Fig. 4.

Fig. B.2 Crosscuts along the planes XY, XZ, and YZ in the 3D dust extinction map of Lallement et al. (2019). We use the Heliocentric Galactic coordinates. The (common) gray scale shows log(δAv), with δAv the differential of the dust extinction in units of magnitude per parsec. The blue and orange dots mark the inner and outer surfaces of the Local Bubble shell measured from the radial profiles of δAv. The purple and brown continuous lines trace the intersection of our models for $r_{\text{inner}}^{{\ell }_{\max }=6}$ and $r_{\text{outer}}^{{\ell }_{\text{max}}=6}$ with the respective plane, respectively.

Fig. B.3 Thickness of the Local Bubble shell as seen from the Sun. The full sky map is shown on the left and latitude profiles for longitudes in steps of 30° are shown on the right.

We point out, however, that this geometric construction is likely to overestimate the distance to the outer surface of the local bubble shell and, consequently, its thickness. This is primarily due to the limited resolution of the 3D dust map and the post-processing of the radial profiles to extract its inflection points. Their effect is that any density peaks are radially smeared out, leading to lower and higher values for the inner and outer radius of the shell. In addition, although the geometric approach is in principle suitable for characterizing the geometry of a shell formed in a homogeneous environment, we have little control over the possible effects of initial inhomogeneities on our results. If, for a given line of sight, the shell merges with an initial overdensity, it is possible that the second inflection point does not correspond to the extreme edge of the shell but to that of the overdensity, or to any intermediate value, depending on the relative amplitudes of the density peaks. In Fig. B.3 we show the thickness of the shell of the Local Bubble as measured from the Sun in the form of a sky map and latitude profiles for several longitudes. It is seen that the thickness varies substantially over the sky and is generally larger in the northern hemisphere than in the southern one. We emphasize the fact that this is the thickness of the shell as measured from the Sun and that, due to projection effects, it is not directly related to the amplification of the magnetic field as it is not measured from the explosion center.

References

All Tables

All Figures

Fig. 1 Crosscuts along the planes XY, XZ, and YZ through a spherical bubble with an inner radius of 100 pc and a thickness of 20 pc. The explosion center (black asterisk) is at the origin of the axes and at the center of the bubble. The direction of B0 is set by $\left({\theta }_{B^0},{\varphi }_{B^0}\right)=\left({60}^{°},0^{°}\right)$ in a spherical coordinate where θ = 0 points to positive Z and (θ, ϕ) = (0, 0) points to positive X. B0 belongs to the XZ plane. The background colors indicate the amplification of the total magnetic field in the plane (|B|/|B0|) and streamline the orientation of the projected B in the crosscut planes. The red arrow in the lower left corner of each panel indicates the projection of B0 in the planes. The red circle has a radius of B0 and indicates the maximum length of the arrow

In the text

	Fig. 2 Dependence of the amplification factor of the magnetic field in the shell of a spherical bubble as a function of the radial distance (top) and as a function of the angle α between the radial unit vector and the direction of the initial magnetic field (bottom), for several values of α and r, respectively. In this example, Rmin = 100 pc and Rmax = 120 pc.
In the text

Fig. 3 Crosscuts along the XZ plane through three bubble cases: a spherical bubble with an off-centered explosion center (left), a nonspherical bubble with the explosion center at the origin of the axes (middle), and that off-centered (right). The explosion center is marked with a black asterisk and the origin of the axes with a black dot. The bubble shape for the middle and right panels are the same. The direction of the initial magnetic field (also marked with the red arrow) is the same for every panel, and as in Fig. 1. B0 belongs to the XZ plane. Colors and streamlines follow the same convention as in Fig. 1.

In the text

Fig. 4 Mollweide projection of the full-sky maps of the contributions from the shell of the Local Bubble to the RM, Q, and U signal (from left to right) as predicted for our different scenarios (rows). Maps are given in Galactic coordinates with the Galactic center at the center of the maps, Galactic longitude increases to the left and the Galactic north pole is at the top. Color scales are shared in columns. They range from -10 to 10 rad/m2 for the RM maps and from –10 to 10 µK for Q and U maps.

In the text

	Fig. 5 Longitude profiles for constant latitude stripes of 10° width for |b| ∈ [55°, 85°], as described in the text. Predictions for the different scenarios for the Local Bubble magnetic field are shown as continuous lines according to the legend, for the RM, Q, and U signal (from left to right). For reference, observational data is also shown as gray points with symmetric errorbars.
In the text

Fig. 6 Orthographic projection of the same data used in Unger & Farrar (2024) (top), the model predictions for the SCA scenario (middle), and pulls (bottom). The gray area is masked out as explained in the text. The north pole is at the center of the left disk and the south pole to the right. The longitude zero is marked by the vertical thick lines starting from the poles upward. Longitude increases clockwise on the left and counterclockwise on the right.

In the text

	Fig. 7 Contribution to the reduced χ2 of the different observables as a function of the strength of the initial magnetic field for the six scenarios. The contribution from Q and U are combined. The contributions from the northern and southern hemispheres are also shown with dashed and dotted lines, respectively.
In the text

	Fig. B.1 Sky maps of the modeled distances to the inner $\left(r_{\text{inner}}^{{\ell }_{\max }=6},\text{left}\right)$ and outer $\left(r_{\text{outer}}^{{\ell }_{\text{max}}=6},\text{right}\right)$ surfaces of the shell of the Local Bubble as measured from the Sun. Projection and conventions are the same as in Fig. 4.
In the text

Fig. B.2 Crosscuts along the planes XY, XZ, and YZ in the 3D dust extinction map of Lallement et al. (2019). We use the Heliocentric Galactic coordinates. The (common) gray scale shows log(δAv), with δAv the differential of the dust extinction in units of magnitude per parsec. The blue and orange dots mark the inner and outer surfaces of the Local Bubble shell measured from the radial profiles of δAv. The purple and brown continuous lines trace the intersection of our models for $r_{\text{inner}}^{{\ell }_{\max }=6}$ and $r_{\text{outer}}^{{\ell }_{\text{max}}=6}$ with the respective plane, respectively.

In the text

	Fig. B.3 Thickness of the Local Bubble shell as seen from the Sun. The full sky map is shown on the left and latitude profiles for longitudes in steps of 30° are shown on the right.
In the text