# Current potential patches

Todd Elder1,2 and Allen H Boozer1 1 Columbia University, New York, NY  2
University of Maryland, College Park, MD

###### Abstract

A novel form of the current potential, a mathematical tool for the design of
stellarators and stellarator coils, is developed. Specifically, these are
current potentials with a finite-element-like basis, called current potential
patches. Current potential patches leverage the relationship between
distributions of magnetic dipoles and current potentials to explore limits of
the access properties of stellarator coil sets. An example calculation is
shown using the Helically Symmetric Experiment (HSX) equilibrium,
demonstrating the method’s use in coil design and understanding the limits of
the access properties of coil sets. Current potential patches have additional
desirable properties such as of promoting sparse current sheet solutions and
identifying crucial locations of shaping current placement. A result is found
for the HSX equilibrium that shaping currents covering only 25% of the winding
surface is sufficient to produce the equilibrium to a good accuracy.

##  I Introduction

Stellarators are devices for magnetic confinement fusion (MCF) which have the
potential to provide nearly limitless clean energy. Modern stellarators are
optimized a priori to their operation using supercomputers and novel
calculational methods, yielding steady-state MCF devices which mitigate
physical problems common to both stellarators and tokamaks such as magnetic
island formation and energetic particle loss. As such, stellarators are at the
cutting edge of MCF technology. The experimental success of the optimized W7-X
stellarator in Germany has confirmed that excellent plasma properties can be
achieved by computational design ConfirmW7XToplogy:2016 .

A major barrier in the development of the stellarator is coil complexity.
Stellarator coils suffer from poor access properties and tight tolerances. In
order for the stellarator to become a device as ubiquitous as the tokamak,
stellarator coils must be simplified. Tokamaks rely on axisymmetry for their
confinement properties. To minimize toroidal ripple, a perfect cladding of
planar toroidal field coils would yield perfect axisymmetry but terrible
access. Alternatively, stellarators can uniquely have far superior access to
the vacuum chamber compared to tokamaks. A method to determine the maximal
access properties for a given equilibrium is developed in this paper.

The access properties of a coil set is typically given by two metrics: coil-
coil spacing and coil-plasma spacing. Coil sets have good access properties,
or are open-access, if the coil-coil spacing is large. Large gaps between
coils allow ease of assembly and access into the vacuum chamber for
maintenance.

In this paper, an alternative formulation of access properties is used: the
fractional area of the winding surface populated by surface currents. A given
equilibrium is found to have good access properties if the fraction of the
winding surface needed to produce the equilibrium is small. The method is a
combination of L0 regularization to find regions of the winding surface where
surface currents are unnecessary and singular value decomposition (SVD) to
ensure the resulting surface currents used are low in amplitude and are placed
in the most efficient locations on the winding surface to support the desired
equilibrium. An example calculation using the HSX equilibrium is shown,
finding that the equilibrium can be supported by placing currents on only 25%
of the HSX winding surface.

###  I.1 Conditions on the coils

Magnetic field coils must satisfy four properties: they must (i) produce
magnetic surfaces (ii) be close enough to the plasma to produce the desired
field (iii) be feasible in an engineering sense and (iv) separated from the
plasma far enough to provide good chamber access. In a fusion power plant,
(iv) is a more important constraint as the coils must be sufficiently
separated from the plasma to allow room for tritium breeding and neutron
shielding.

####  I.1.1 Formation of nested magnetic surfaces

The coil currents and geometries must produce a magnetic surface with
$\mathbf{B}\cdot\hat{n}=0$ on the surface.

Magnetic confinement uses toroidal magnetic surfaces. Toroidal magnetic
surfaces are generated by magnetic field lines which wind around the torus
without leaving the specified surface, $\mathbf{B}\cdot\hat{n}=0$ where
$\hat{n}$ is the unit normal to the surface. The magnetic field has
contributions both from external currents and the plasma current of the
device.

For axisymmetric devices, these surfaces are crudely in the shape of a donut.
For non-axisymmetric devices like stellarators, the surfaces have a
complicated three-dimensional geometry which varies along the toroidal
direction.

####  I.1.2 Close enough to the surface to produce the required fields

Away from the coils the magnetic field is curl-free, $\mathbf{B}=\nabla\Phi$.
This, combined with the divergence-free condition of magnetic fields yields
Laplace’s equation $\nabla^{2}\Phi=0$ which in a cylinder which is periodic
along the $\hat{z}$ axis of the cylinder - an approximate geometry of a torus
- has solutions which depend exponentially on the radial coordinate. This
simple fact defines limits for magnetic surface shaping, as the spatial decay
of magnetic fields produced by distant currents defines the limits of feasible
surface shaping for stellarators. In particular, only magnetic fields which
decay slowly through space when produced by a current on the cylinder may be
used to design stellarator reactors. These fields are known as the efficient
field distributions Boozer:Stell_Design  EFDs:2016 .

####  I.1.3 Feasibility

The coils cannot have (i) current densities that are beyond what is producible
by coils (ii) forces on them that are too large (iii) curvatures or torsions
that are too great and (iv) must be sufficiently widely separated.

####  I.1.4 Provide adequate plasma chamber access

It is important to have good chamber access for ease of assembly and
maintenance. For example, the small plasma-chamber distance of W7-X was a
major contributor to budget overruns and schedule delays ZhuSimplerCoils .
Designing coils that optimize access is a focus of this work. Access is
quantified by the fraction of the gridded winding surface used by shaping
currents of current potential patches - the less of the winding surface used,
the greater the access properties.

##  II Background

The initial configuration space for 3D coils to produce a given magnetic
surface is too large to search by brute force. Good coil guesses can be found
using physics motivations. A typical method of doing this was invented by
Rehker and Wobig ModCoils_1st:1973  in which a current sheet near the desired
magnetic surface is found.

The formalism of these magnetic-surface-generating current sheets was
developed by Merkel NESCOIL:1987  and used to successfully develop coil sets
for W7-AS and HSX. Merkel extended the winding surface from infinitesimally
separated, as in the Rehker and Wobig formulation, to a surface placed an
arbitrary distance away from the desired magnetic surface. Merkel posed the
problem as solving for what is now known as the current potential $\Phi$.
Details of the formulation may be found In Landreman 2017 REGCOIL:2017 .

$\Phi$ is known as the current potential, the magnetic surface generating
stream function.

$\Phi(\theta,\zeta)=I_{\theta}\zeta/2\pi+I_{\zeta}\theta/2\pi+\Phi_{SV}(\theta,\zeta)$
(1)

and has two terms: the non-secular and secular (also known as single-valued
and multi-valued) components. The secular components are
$I_{\theta}(\zeta/2\pi)+I_{\zeta}(\theta/2\pi)$ which increase unbounded with
poloidal and toroidal transits by $2\pi$. The non-secular components

$\Phi_{SV}(\theta,\zeta)=\sum_{m,n}\Phi_{m,n}^{c}\cos(m\theta-n\zeta)+\Phi_{m,n}^{s}\sin(m\theta-n\zeta)$
(2)

are single-valued functions, typically in the Fourier basis. The non-secular
terms are degrees of freedom for generating the desired magnetic field. The
secular terms are degrees of freedom for the current potential. $I_{p}$ is
fixed by the toroidal flux passing through the plasma. $I_{t}$ is a degree of
freedom to control whether the current potential is helical or not, and
provides the magnetic flux passing through the hole in the torus. If $I_{p}=0$
the current potential is unable to produce toroidal flux. If $I_{p},I_{t}=0$
the current potential is single-valued.

A prescribed magnetic surface is generated with currents by eliminating normal
components of the magnetic field $\mathbf{B}$ that arise due to its
interaction with the surface geometry. $\mathbf{B}$ can have many sources,
such as

$\mathbf{B}=\mathbf{B}_{\mbox{plasma}}+\mathbf{B}_{\mbox{external}}$

where $\mathbf{B}_{\mbox{plasma}}$ is the plasma’s contribution to
$\mathbf{B}$, yielding a solenoidal field, and $\mathbf{B}_{\mbox{external}}$
is the field generated by an external coil set.

For example, if designing a coil set consisting of windowpane coils and
toroidal field coils, one would begin by computing the
$\mathbf{B}_{n}=\mathbf{B}_{\mbox{coils}}\cdot\hat{n}$ and use this to compute
a single-valued current potential to generate windowpane contours.

Solving for $\Phi_{SV}$ is accomplished using an integral formulation of
$\mathbf{B}\cdot\hat{n}=0$, see Landreman 2017 REGCOIL:2017  for more details.
For the purposes of this work we restrict ourselves to only considering
sources of $\mathbf{B}\cdot\hat{n}$ from vacuum magnetic fields
$\mathbf{B}=\nabla\Phi$ which is a purely toroidal $1/R$ field. The problem
then becomes a typical least-squares problem to minimize

$\displaystyle\min\bigg{(}\big{|}\mathbf{A}\cdot\mathbf{x}-\mathbf{b}\big{|}_{2}\bigg{)}$
(3)

where $\mathbf{A}$ is known as the transfer or inductance matrix, the
components of $\mathbf{b}$ are the values of $\mathbf{B}\cdot\hat{n}$ for each
point on the surface $S$, and the components of $\mathbf{x}$ are the degrees
of freedom of the current potential to satisfy $\mathbf{B}\cdot\hat{n}=0$ on
the surface. In the Fourier basis, the components of $\mathbf{x}$ correspond
to the amplitudes $\Phi_{m,n}^{c,s}$.

The distinction between transfer and inductance matrices lies in the source of
the field EFDs:2016  traveling from the winding surface to the magnetic
surface. Transfer matrices describe the spatial decay of a field from the
winding surface to the magnetic surface, for example due to a point dipole
source. Inductance matrices relate currents to the field they produce. For
current potentials, $\mathbf{A}$ relates the current on the winding surface to
the magnetic field produced on the surface $S$. $\mathbf{A}$ relates the
geometry of the two surfaces involved, $S$ and $S^{\prime}$, and the gradient
of the current potential $\nabla\Phi$. The inductance matrix can be very ill-
conditioned and so represents the main computational difficulty of this least-
squares problem. More on the inductance matrix can be found here Lazanja:2011
. Further insights and a rapid way to compute the inductance matrix may be
found here Cerfon_Multipole:2012 .

####  II.0.1 Relationship to distributions of magnetic dipoles

An important but underappreciated result of Merkel’s work is the relationship
between the current potential and distributions of magnetic dipoles, which we
utilize in this work. This result has been more rigorously derived, for
example in Zhu_PMs:2020  and an alternative derivation included in appendix
Appendix:1.

Equation (6) of Merkel’s work NESCOIL:1987  reads,

$g(\theta^{\prime},\zeta^{\prime},\theta,\zeta)=-\dfrac{1}{4\pi}\sum_{l=0}^{n_{p}-1}\bigg{[}\dfrac{\mathbf{n}\cdot\mathbf{n}^{\prime}}{|\mathbf{x}-\mathbf{x}^{\prime}|^{3}}-\dfrac{3\mathbf{n}\cdot(\mathbf{x}-\mathbf{x}^{\prime})\mathbf{n}^{\prime}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}{|\mathbf{x}-\mathbf{x}^{\prime}|^{5}}\bigg{]}$
(4)

where $g$ is the weight for each Fourier harmonic of the current potential in
the inductance matrix, $\mathbf{x}(\theta,\zeta)$,
$\mathbf{x}^{\prime}(\theta^{\prime},\zeta^{\prime})$ are the position vectors
on the winding and plasma surfaces, and $\mathbf{x}(\theta,\zeta)$,
$\mathbf{x}^{\prime}(\theta^{\prime},\zeta^{\prime})$ are the unit normal
vectors. Compared to the field of a magnetic dipole,

$\mathbf{B}(\mathbf{x})=\dfrac{\mu_{0}}{4\pi}\dfrac{3(\mathbf{x}-\mathbf{x}^{\prime})((\mathbf{x}-\mathbf{x}^{\prime})\cdot\mathbf{m})-\mathbf{m}}{|\mathbf{x}-\mathbf{x}^{\prime}|^{3}}$
(5)

where $\mathbf{m}$ is the magnetic moment, the relation between the two
formulations is apparent. This relationship has been expanded upon by Boozer
Boozer:Stell_Design EFDs:2016  in the efficient field distributions framework.
This relationship between magnetic dipoles and current potentials also
provides a simple basis for creating discontinuous patches of current sheets.
A simple derivation is given in appendix Appendix:1.

##  III Method

Figure 1: Example configuration of gridded magnetic dipoles on the HSX winding
surface. Each dipole amplitude has been optimized to create the HSX magnetic
surface. Only one field period of both the winding and magnetic surface is
necessary to capture the full toroidal problem. For this case, the magnetic
surface is produced with a maximal field error of $10^{-5}$%. The total field
contributing to $B_{n}$ here is the toroidal $1/R$ field and the field of the
magnetic dipoles.

Magnetic dipoles are placed in a grid on the HSX winding surface (see Fig.
fig:3D_Dipoles+MagSurf) and are oriented normal to the winding surface. The
winding surface used is the surface on which HSX’s modular coils lie. The
magnetic surface used is a truncated version of HSX’s experimental
equilibrium. The truncation minimizes the effect of modular coil ripple
present in the equilibrium arising from the 12 coils per period of HSX.
However, as discussed in Fig. fig:Bn_Fourier, this method does not entirely
eliminate defects from coil ripple.

Each dipole acts as a local element to contribute magnetic field. The dipole
strengths are then adjusted via a least-squares problem to minimize the normal
field on the surface due to the a background toroidal field,

$\min\bigg{(}|\mathbf{A}_{\mbox{dip}}\cdot\mathbf{x}-\mathbf{b}|^{2}\bigg{)}$
(6)

where $\mathbf{A}_{\mbox{dip}}$, defined below, is the transfer matrix of the
dipole magnetic fields, $\mathbf{x}$ represents the individual dipole
strengths, and $\mathbf{b}$ represents the normal magnetic field at each point
on the magnetic surface produced by external coils. In later steps, dipoles
are removed entirely from the winding surface to enforce solution sparsity, or
the open-access properties of the coil set.

Other sources of $B_{n}$ than a purely toroidal field were considered, for
example from simple helical coils, though the results did not change
significantly, neither qualitatively or numerically, and so only results from
$B_{n}$ arising from a purely toroidal field are shown.

The method is motivated by the idea of a current potential with regions of
high current density in some locations and low in others, or a ”sparse current
potential.” Sparse current potentials are intended to have concentrated
currents in the sheet which lie in small regions or thin bands. These features
cannot easily be realized in a Fourier basis due to its global nature. A
finite element method (FEM) basis would be ideal, in which localized patches
of the current sheet can contribute independently. FEM bases support feature
sharpness and so are amenable to sparsity.

Magnetic dipoles are theoretically rings of current of vanishing radius. As a
result, the current of this basis is trivially divergence-free, however
current cannot flow between patches. Indeed, each dipole is not a patch of
current, but rather a $\delta$ function of current centered in each patch.
This is the approximate nature of this basis. The surface current flowing on
the winding surface is given exactly by the gradient of the current potential,
which may be calculated using standard numerical methods over a grid. This
approach is similar to the work done by Helander Helander_PMs:2020 . Despite
this approximation, the magnetic dipole basis current potentials found agree
well with Fourier basis current potentials for the same magnetic surfaces, see
Sec.sec:REGCOIL_comparison.

The use of magnetic dipoles, specifically permanent magnets, to generate
shaping fields is a recent approach to stellarator design Helander_PMs:2020 .
In this approach, a surface or volume grid of magnetic dipoles is used to
generate a given magnetic surface by varying dipole strength and orientation.
These magnetic dipole methods have been improved and implemented for
optimization Zhu_PMs:2020  Kaptanoglu_GPMO:2023  Hammond_PMs:2020  leading to
the construction of the first permanent-magnet stellarator MUSE MUSE:2021 .
Permanent magnet stellarators are remarkably able to produce magnetic surfaces
to great accuracy and provide some open access, though are not viable for
fusion power plants due to a variety of reasons: degradation of the magnets
from neutron flux and the low field strength of the magnets if they are placed
behind shielding.

Here we will describe how to initialize the magnetic dipole configurations
using the winding surface. The winding surface may be represented as a surface
$\mathbf{S}_{w}(\theta,\zeta)$ with a mesh grid separated uniformly in
$\theta$ and $\zeta$. The magnetic dipole winding surface then consists of
$n_{\theta}$ by $n_{\zeta}$ points, each of which is populated by a single
magnetic dipole.

The field generated by each magnetic dipole on a plasma magnetic surface
$\mathbf{S}_{p}$ at the $i$th point on the surface $\mathbf{S}_{p,i}$ is

$\mathbf{B}(\mathbf{S}_{p,i})=\dfrac{3\hat{\mathbf{r}}(\mathbf{m}_{j}\cdot\hat{\mathbf{r}}-\mathbf{m}_{j})}{|\mathbf{r}|^{3}}$

where $\mathbf{r}=\mathbf{x}_{j}-\mathbf{S}_{p,i}$,
$\hat{\mathbf{r}}=\mathbf{r}/|\mathbf{r}|$, $\mathbf{x}_{j}$ is the position
of the $j$th dipole, and $\mathbf{m}_{j}$ is the $j$th magnetic dipole vector.
The inductance matrix is accordingly

$\mathbf{A}_{\mbox{dip}}=\begin{bmatrix}\dfrac{3\hat{\mathbf{r}}_{11}(\mathbf{m}_{1}\cdot\hat{\mathbf{r}}_{11}-\mathbf{m}_{1})}{|\mathbf{r}_{11}|^{3}}&\dfrac{3\hat{\mathbf{r}}_{21}(\mathbf{m}_{2}\cdot\hat{\mathbf{r}}_{21}-\mathbf{m}_{2})}{|\mathbf{r}_{21}|^{3}}&...\\\
\dfrac{3\hat{\mathbf{r}}_{12}(\mathbf{m}_{1}\cdot\hat{\mathbf{r}}_{12}-\mathbf{m}_{1})}{|\mathbf{r}_{12}|^{3}}&...&...\\\
...&...&...\\\ \end{bmatrix}$

which is an $N\times M$ matrix where $N$ is the number of evaluation points on
the magnetic surface and $M$ is the number of dipoles used. Stellarator
symmetry may be used to account for each dipole’s contribution from being
mapped toroidally. As a result, only a half-field-period winding surface with
a half-field-period of plasma surface need be considered when solving for the
magnetic dipole strengths.

##  IV Results

Figure 2: Surface plot of magnetic dipole amplitudes $m$ on one half field
period of the HSX winding surface for a background toroidal field, analogous
to the current potential $\Phi$. From left to right the dipole amplitudes are
SVD-regularized with cutoff threshold $10^{-5}$, $10^{-10}$, $10^{-15}$
corresponding to limits on the maximum dipole amplitudes. For reference, SVD-
regularization of the order $10^{-5}$ corresponds to dipoles 10x stronger than
permanent magnets. pRes refers to plasma surface resolution, while wRes refers
to winding surface resolution. From top to bottom the dipole density increases
as grids of 128x64, 160x80, and 192x86 magnetic dipoles. The highest-
resolution cases best-represent an actual current potential. In the bottom
three figures, we see high-current-density regions appearing as large
windowpane contours. Current flows along contours of constant current
potential so provide initial guesses for windowpane coils.

###  IV.1 Winding surface solutions without sparsity promotion

The density of magnetic dipoles on the surface influences the error field
significantly. As we populate the winding surface with more and more dipoles
($n\to\infty$), the discrete set of dipoles approximates a continuous FEM
basis. The effect of this increasing density can also be seen in Fig.
fig:Dipole_CPs_Full_WS. The physical interpretation of this effect is that
higher dipole densities can support shorter-wavelength Fourier spectra.

To begin, eqn. eqn:CP_LS_problem is solved using the full HSX winding surface
with different dipole densities on the surface $n_{d}$ and levels of singular
value decomposition (SVD) truncation. Here, SVD regularization is the
suppression of singular values of eqn:CP_LS_problem less than some cutoff
parameter $r_{\mbox{cond}}$. The effect of these two parameters may be seen in
Fig. fig:Dipole_CPs_Full_WS. Greater dipole densities $n_{d}$ allow the
formation of high-wavenumber features in the current potential. SVD
regularization limits the maximal amplitude of individual dipoles and forces
solutions to be smooth by cutting off inefficient modes, see Sec.
sec:SVD_fractional_regularization for more on this.

A closer inspection of the condition number of the inductance matrix reveals
the tendency for high-m harmonics to arise in the current potential. The
condition number of eqn:Dipole_inductance_matrix when using a grid of 128
toroidal cells by 192 poloidal cells is $\sim 10^{20}$, suggesting these
solutions are very sensitive to perturbations and so the problem must be posed
in such a way that the condition number is small. Thus, while these full-
winding-surface-solves reveal qualities of the dipole basis, feasible coil
sets must be generated using a direct SVD-truncation technique, Sec.
sec:SVD_fractional_regularization.

####  IV.1.1 Sparsity through $L0$ regularization and SVD truncation of
current potential patches

Figure 3: Fourier decomposition of Bn on the magnetic surface due to a
toroidal flux. At 0.1% field accuracy, the magnetic surface quality becomes
”good enough” as a rough measurement. Interestingly, the HSX equilibrium used
has some very high mode number harmonics, likely due to this equilibrium not
being the original stage-1 result of HSX which has been lost. The importance
of this high-mode harmonic to the equilibrium recreation could introduce
significant complexity into the corresponding current configurations. Figure
4: Condition number of the dipole-basis transfer matrix. As expected, it
increases exponentially with the number of singular values included. Condition
numbers between 5-10 are expected to give feasible current distributions which
produce the field with sufficient accuracy to recover solution sparsity.
Figure 5: SVD-truncated dipole-basis current potentials. The chosen condition
numbers are 10,20, and 30. Solution sparsity is achievable but there is
interestingly still high ripple in the solutions. However, these features
correspond to low current density regions and so are perhaps more feasible
current arrangements. Conversely, current potentials with high condition
numbers involve large regions of high current density, a less feasible current
arrangement. Figure 6: Fraction of surface area required to produce the
equilibrium to 0.1% accuracy for different condition numbers of M. As
complexity of M increases, more solution sparsity emerges. However, the
saturation of the fractional surface required to produce good-quality magnetic
surfaces with respect to matrix condition number provides estimate limits for
open-access properties of the coil set. The level of access achieved before
saturation is unique to each equilibrium.

The method that follows is closely related to the efficient field
distributions method of EFDs:2016 , though this implementation uses the unique
ability of current potential patches to construct discontinuous current
sheets. A given normal magnetic field $B_{n}(\mathbf{x})$ on the plasma
surface may be Fourier decomposed as $B_{n}(\mathbf{x})=b_{k}(\theta,\zeta)$
where

$b_{k}(\theta,\zeta)=\sum_{m,n}^{k}b_{m,n}\cos(m\theta-n\zeta)$ (7)

where the truncation value $k$ is determined by the maximum tolerable field
error.

Exponential convergence of the Fourier series allows us to choose a Fourier
truncation value $k$ such that the maximal error $\delta$ between the normal
field and its Fourier representation is $|B_{n}-b_{k}|^{2}<\delta$. This
Fourier representation for a certain $\delta$-quality magnetic surface tells
us the efficiency with which it can be produced by a grid of magnetic dipoles
on some external surface,

$\displaystyle\mathbf{M}\cdot\mathbf{d}=\mathbf{b}_{k}$ (8)

and the singular value decomposition of $\mathbf{M}$ informs the efficiency of
the dipoles,

$\displaystyle\mathbf{M}=\mathbf{U}\cdot\mathbf{S}\cdot\mathbf{V}^{\dagger}$
(9)

$\mathbf{U}^{\dagger}\mathbf{d}$ relates which distributions of dipoles are
most efficient at producing Fourier components, while
$\mathbf{V}^{\dagger}\cdot\mathbf{b}_{k}$ relates which magnetic fourier terms
are easily produced by the dipoles. $\mathbf{M}$ is similar to the transfer
matrix $\mathbf{A}$ but instead of relating the magnetic field of magnetic
dipoles to points on the magnetic surface, $\mathbf{M}$ relates point dipoles
to Fourier-decomposed normal magnetic field distributions on the magnetic
surface.

The condition number of $\mathbf{M}$ is given by the ratio of the largest to
smallest singular values of $\mathbf{S}$. If the condition number is high,
then inefficient components of the dipole field are being used to produce the
field on the surface. As long as the condition number is low, say
$10^{0}-10^{1}$, the dipoles are used efficiently and the resulting current
potentials are more feasible.

We apply this analysis technique to the HSX equilibrium. The Fourier
representation of $B_{n}$ eqn:Bn_Fourier may be seen in Fig. fig:Bn_Fourier.
The condition number of the associated matrix $\mathbf{M}$ from eqn:M_SVD may
be found in Fig. fig:CondNumber.

We enforce sparsity through the following method: eqn:M_SVD is truncated of
singular values up to a given condition number, which is then solved via a
least-squares problem eqn:Dipole_Bn_efficiency_discussion, then sparsity is
enforced by removing all dipoles with strengths less than some threshold
strength $m_{dip}$ from the winding surface. Finally, the remaining dipole
strengths are adjusted through another least-squares solved using only the
remaining dipoles place on the winding surface. In this way, we may use small
coverage regions of the current potential surface to produce the desired
equilibrium while maintaining a low matrix condition number.

Example truncated-SVD current potentials may be found in Fig. fig:Money1. The
more efficient locations of current placement coincide wihh low current
density regions of the current potential. In Fig. fig:Dipole_CPs_Full_WS,
these are the whitish regions of the current potential. The implication
coincides with a general feature of efficient fields: inefficient fields are
produced using short-wavelength, high current density features in the current
potential.

###  IV.2 Comparison with REGCOIL

We compare the current potential patch results with single-valued current
potentials in the Fourier basis, fig:REGCOIL_SV_CPs. Each current potential
shown produces the HSX magnetic surface to similar levels of error to the
figures found in Fig. fig:Dipole_CPs_Full_WS using different levels of
Tikhonov regularization. In this case, Tikhonov regularization uses the $L2$
norm of the current density, which is the REGCOIL formulation REGCOIL:2017 .
The Fourier-basis current potentials do remarkably well in producing sparse
features, likely due to the high toroidal and poloidal mode numbers allowed in
the current potential. Even so, these results suggest strong localization of
current can be achieved with fourier-basis current potentials.

##  V Summary

A novel method to investigate access properties of stellarator coil sets was
developed. The tool, known as current potential patches, uniquely retain the
solution properties of current potentials for coil design yet are able to
generate sparse solutions. An example calculation using the HSX winding and
magnetic surfaces demonstrates the ability of current potential patches to (i)
identify crucial locations of current placement for surface shaping, for
example using windowpane coils and (ii) provide an example bounding
calculation on a coil set with maximal open-access properties which still
produces the desired magnetic surface to a given degree of accuracy. This
method could be applied to (iii) the generation of localized windowpane coils
for error field remediation or the use of a coilset consisting of simple
modular or helical coils to provide the bulk toroidal field and windowpane
coils to provide the rest. Overall, this method has use in stellarator design
in which access to the vacuum vessel and region interior to the coil set is
desired.

The access limits set by the method are approximate due to the use of surface
currents. The use of filamentary coils would likely decrease access properties
due to the inability of filamentary coils to produce the current distributions
of current sheets.

This method has a few direct applications. For one, coil sets consisting of
simple toroidal field coils or helical coils in conjunction with windowpane
coils can be generated using this method. In this case, contours of
$\Phi_{SV}$ are the contours of windowpane coils. Second, a metric describing
the amenability of an equilibrium to open-access coil sets could be derived.
Further, equilibria which support good access properties could be simpler by
not relying on complicated shaping currents throughout the winding surface to
support them.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Fusion Energy Sciences under Awards DE-
FG02-95ER54333 and DE-FG02-03ER54696.

##  Appendix A Single-valued current potential and dipoles

The single-valued part of the current potential $\Phi$ in a small region on a
surface can be represented by a magnetic dipole. Let $\hat{z}$ be the normal
to the surface and use cylindrical coordinates. The magnetic moment analog of
a current density is

$\displaystyle\vec{m}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
r\hat{r}\times\vec{j}d^{3}x,\mbox{ and }$ (10) $\displaystyle\vec{j}$
$\displaystyle=$
$\displaystyle(\vec{\nabla}\kappa\times\hat{z})\delta(z),\mbox{ so }$ (11)
$\displaystyle\vec{m}$ $\displaystyle=$ $\displaystyle-\frac{\hat{z}}{2}\int
r\frac{\partial\kappa(r,\theta)}{\partial r}\delta(z)rdrd\theta dz$ (12)
$\displaystyle=$ $\displaystyle-\frac{\hat{z}}{2}\int
r^{2}\left(\int_{0}^{2\pi}\frac{\partial\kappa(r,\theta)}{\partial
r}d\theta\right)dr$ (13) $\displaystyle=$ $\displaystyle-\frac{\hat{z}}{2}\int
2\pi r^{2}\bar{\kappa}(r)dr$ (14) $\displaystyle=$
$\displaystyle-\pi\hat{z}\int\Big{(}\frac{d}{dr}(r^{2}\bar{\kappa})-2r\bar{\kappa}dr$
(15) $\displaystyle=$ $\displaystyle\hat{z}\int\kappa(r,\theta)rdrd\theta,$
(16)

where $2\pi\bar{\kappa}\equiv\int_{0}^{2\pi}\kappa(r,\theta)d\theta$. Since
$(\hat{r}\times\hat{\theta})\cdot\hat{z}=1$,
$\hat{r}\times(\vec{\nabla}\kappa\times\hat{z})=-\hat{z}\partial\kappa/\partial
r$. The current potential $\kappa(r,\theta)$ is assumed to be zero outside of
the region of the radial integration. The area element for this patch on the
surface is $da=rdrd\theta$.

Figure 7: Single-valued fourier-basis current potentials computed by REGCOIL,
each with a different level of Tikhonov regularization. The top figures are
$\Phi_{SV}$ while the bottom are $|K|$. The title of each figure denotes the
integrated current density $|K|$, amount of Tikhonov regularization $\lambda$,
and average $|B_{n}|$ error. These current potentials were generated using a
grid of 64 poloidal points and 32 toroidal points per half field period, and a
maximum poloidal and toroidal mode number of 14. Remarkably, the dipole basis
produces similar qualitative results to that of Fig. 2: the high current
density locations are similarly localized and have high ripple.

## References

  * (1) S. I. Krasheninnikov and A. S. Kukushkin, _Physics of ultimate detachment of a tokamak divertor plasma_ , Journal of Plasma Physics, 83, 155830501 (2017).
  * (2) TS Pedersen et. al., _Confirmation of the topology of the Wendelstein 7-X magnetic field to better than 1:100,000_ , Nature communications, 7 1, 13493 (2016).
  * (3) C Zhu et. al., _Towards simpler coils for optimized stellarators_ , IAEA Conference, (2018).
  * (4) S. Rehker and H. Wobig, _Stellarator Fields with Twisted Coils_ , IPP 2/215 , (1973).
  * (5) P. Merkel, _Solution of stellarator boundary value problems with external currents_ , Nuclear Fusion, 27, 867 (1987).
  * (6) M. Landreman, _An improved current potential method for fast computation of stellarator coil shapes_ , Nuclear Fusion, 57, 046003 (2017).
  * (7) Allen H Boozer, _Physics of magnetically confined plasmas_ , Reviews of modern physics, 76 4, 1071-1141 (2005)
  * (8) AH Boozer, _Stellarator design_ , Journal of Plasma Physics, 81 6, 515810606 (2015).
  * (9) M. Landreman and AH Boozer, _Efficient magnetic fields for supporting toroidal plasmas_ , Physics of Plasmas, 23, 032506 (2016).
  * (10) D. Lazanja, _Components of the magnetic field from interior and exterior sources in toroidal geometry_ , Columbia University PhD thesis, (2011).
  * (11) T. Askham and A.J.Cerfon, _An adaptive fast multipole accelerated Poisson solver for complex geometries_ , Journal of Computational Physics, 304, 1-22 (2017).
  * (12) Caoxiang Zhu et al Designing stellarators using perpendicular permanent magnets Nucl. Fusion 60, 076016 (2020).
  * (13) P. Helander et. al., _Stellarators with Permanent Magnets_ , Physical Review Letters, 124, 095001 (2020).
  * (14) AA Kaptanoglu, GP Langlois, M Landreman, _Topology optimization for inverse magnetostatics as sparse regression: application to electromagnetic coils for stellarators_ , Computer Methods in Applied Mechanics and Engineering, 418A, 116504 (2023).
  * (15) K.C. Hammond et. al., _Geometric concepts for stellarator permanent magnet arrays_ , Nuclear Fusion, 60, 106010 (2020).
  * (16) M. Zarnstorff et. al., _MUSE: A Simple Optimized Stellarator Using Permanent Magnets_ , APS DPP conference (2021).


# Interaction of the Galactic-Centre Super Bubbles with the Gaseous Disc

Yoshiaki Sofue1 and Jun Kataoka2
1\. Institute of Astronomy, The University of Tokyo, Mitaka, Tokyo 181-0015,
Japan
2\. Faculty of Science and Engineering, Waseda University, Shinjyuku, Tokyo,
169-8555, Japan E-mail<EMAIL_ADDRESS>
###### Abstract

The interaction of Galactic-Centre (GC) super bubbles (GSB) with the gaseous
disc and halo of the Milky Way is investigated using radio continuum, X-ray,
HI and CO line surveys. The radio North Polar Spur (NPS) constitutes the
brightest eastern ridge of GSB, brightening towards the galactic plane and
reaching $l=22^{\circ},\ b=+2^{\circ}$ at the sharpest end, where it
intersects the tangential direction of the 3-kpc expanding ring and crater.
Examination of the spur ridges reveals that the entire GSB, including the NPS
and its counter spurs, constitutes a GC-symmetrical $\Omega/$ $\Omega$ shape.
The thickness and gas density of the HI and CO discs are shown to increase
sharply from the inside (lower longitude) to the outside of the 3-kpc crater.
Formation of crater is explained by the sweeping of the upper layer of disc
gas by the shock wave from the GC by the explosion $\sim 10$ My ago with the
emitted energy of several $10^{55}$ ergs. Based on the discussion, a unified
view on the structure and formation mechanism of GSB is presented.

###### keywords:

ISM: individual objects: (North Polar Spur) – ISM: shock wave – ISM: bubbles –
Galaxy: centre – galaxies: individual: objects (the Milky Way)

##  1 Introduction

The Galactic Centre (GC) exhibits a variety of explosive phenomena, associated
with high-energy outflows and shock waves. The most energetic and gigantic
shock wave so far reported is the GC super bubbles (GSB)111The same object has
been given various names such as north/south polar spurs (NPS/SPS); radio
loops; Galactic super bubbles; bipolar hyper shells (BHS); ROSAT/eROSITA
bubbles, etc..  due to a nuclear explosion with released energy on the order
of $\sim 10^{55-56}$ erg and extent as large as $\sim\pm 10$ kpc above and
below the Galactic plane (Sofue1977, 2000; Sofue et al 2016; Bland-Hauthorn et
al. 2003; Crocker et al. 2015; Sakar et al. 2016; Kataoka et al. 2018). In
this paper, we revisit this prominent object in the Galaxy in view of the GC
explosion model, being inspired by the recent excellent X-ray observations of
the whole sky with the eROSITA [Predehl et al. 2020].

High-altitude, large scale bubbles in the halo around the GC have been
evidenced in the last decades by the all sky surveys in radio continuum
emission [Haslam et al. 1982, Reich et al. 2001, Planck collaboration 2013],
X-rays [Snowden et al. 1997, Predehl et al. 2020], and $\gamma$-rays [Su et
al. 2010], although their scale size and morphology differ in different energy
bands. Moreover, giant magnetized outflows from the GC are also suggested by
linear polarization measurement in radio band [Carretti et al. 2013].

The bipolar double horn structure has been interpreted as a shock wave from
the GC, and is well simulated by the bipolar-hyper shell (BHS) model [Sofue
1994, Sofue 2000, Guo & Mathews 2012, Sofue et al. 2016, Sarkar et al. 2016,
Sarkar 2019]. The connection of the GSB to the energetic phenomena in the GC
has been firmly confirmed by the discovery of the Fermi bubbles [Su et al.
2010, Crocker et al. 2015, Kataoka et al. 2018]. In fact, various outflow
phenomena have been known in the GC region, such as the GC lobe (GCL) [Sofue &
Handa 1984] as recently revealed of its vertical connection to higher latitude
halo [Heywood et al.2019]. An X-ray chimney extending hundreds of parsecs
above and below the GC has been discussed with tight correlation with the
radio and infrared (IR) vertical structures [Nakashima et al. 2013, Ponti et
al. 2019, Ponti et al. 2021]. Disturbed ionized gas appears to be flowing out
at high speed from the GC into the halo below the galactic plane [Savage et
al. 2017]. These observations seem to connect Sgr A to the Fermi bubbles,
further corroborating the GC-origin scenario of the outer bubbles.
Interconnections and coherent features in radio, IR, X- and $\gamma$-ray maps
are thus suggested in the central few degrees of the GC to $\sim 10$ kpc of
the Galactic halo in the past $\sim 10-20$ My.

Fig. 1 summarizes the starburst or GC explosion model [Sofue 2000, Sofue et
al. 2016], which considered the shock wave propagation in the disc and halo,
as compared with the X-ray all-sky maps.

The top panel shows comparison of the model with the ROSAT all-sky map
[Snowden et al. 1997], where the interstellar absorption is calculated for a
disc model of neutral gas. The middle panel shows a simulated three-color
X-ray image in the 0.25, 0.75 and 2 keV compared with the eROSITA observation
for the inner region of GSB [Predehl et al. 2020]. The interstellar absorption
has been corrected for using the HI survey data [Kalberla et al. 2005]. The
simulations well reproduce the observed X-ray morphology of GSB and the
increasing absorption and hardening of the spectrum toward the Galactic plane.

Figure 1:  Summary of the bipolar hyper-shell model. [Top] Calculated X-ray
sky (Sofue 2000) compared with the ROSAT all sky X-ray map (Ssnowden et al.
1997). RGB colors stand for 1/4, 3/4, and 1.5 keV, respectively. For
intensities see the cited papers. [Middle] Simulated X-ray map on the
$100^{\circ}\times 100^{\circ}$ sky (RGB for 0.25, 0.75, and 2 keV in
arbitrary intensity scaling) corrected for extinction (Sofue et al. 2016)
using the HI survey data (Kalberla et al. 2005). [Bottom] eROSITA RGB map of
the central $100^{\circ}\times 100^{\circ}$ region (Predehl et al. 2020) as
reproduced from url https://www.mpe.mpg.de/7461761/news20200619. Note the
significant absorption near the galactic plane, where the radio data are
crucial in order to see the interaction of the spurs with the disc.

Besides the bubbles and outflows, a giant HI hole has been known around the
GC, which composes a bipolar conical cavity of neutral gas and is supposed to
have been blown off by a galactic wind [Lockman 1984, Lockman & McClure-
Griffiths 2016]. The HI hole has its root in the Galactic plane, coinciding
with the 3-kpc expanding ring [Sanders & Prendergast 1974, Bania 1980, Cohen &
Davies 1976]. The cylindrical boundary of the HI hole makes a crater-like
structure of radius 3 kpc around the GC, whose tangential directions coincides
with the roots of the North Polar Spur and its counter spurs [Sofue 2017].

The North Polar Spur (NPS) is the most prominent spur in the galactic halo in
radio [Haslam et al. 1982] and soft X rays [Snowden et al. 1997, Predehl et
al. 2020]. The radio NPS emerges from the galactic plane at G22+02 [Sofue &
Reich 1979], where the spur exhibits maximum brightness and the sharpest cross
section. This position exactly coincides with the tangential direction of the
3-kpc expanding ring in HI and CO lines [Sofue 2017].

Energetic events in the GC associated with giant bubbles and shock waves,
thus, exhibit various foot prints in the halo and disk from neutral ISM to
high-energy emissions. In this paper we re-examine the morphological
correlation between the 3-kpc crater and the GSB, looking for further evidence
of the GC explosion in the HI and molecular discs using radio continuum,
X-ray, HI and CO line surveys. We, then, attempt to give a comprehensive
explanation for various diverse features observed over the wide range of
spectrum.

##  2 X-Ray GSB absorbed by the disc

###  2.1 Extinction of X-ray spurs by the galactic disc

The top panels of Fig. 1 demonstrate the heavy extinction of the X-rays from
the GC and halo by the galactic disc. Fig. 2 shows the latitudinal variation
of the X-ray intensity at 0.75 keV as observed by ROSAT [Snowden et al. 1997].
The X-ray intensities along the spurs are significantly depressed at latitudes
lower than $|b|\lesssim 10^{\circ}$. Such depression can be understood as due
to the following two kinds of bounded features at low latitudes.

(a) Absorption-bounded cut due to interstellar extinction by heavy elements,
whose optical depth is proportional to the line of sight depth of the disc
obeying the sec $|b|$ law [Sofue 1994, Sofue 2015, Lallement et al. 2016]
(Fig. 2). This cut occurs at $|b|<\sim 10-20^{\circ}$ by the local HI disc.

(b) Emission-bounded cut due to suppression of the intrinsic emissivity at low
latitudes according to decreasing temperature in the galactic disc, where
neutral gas is dominant at temperatures below $<\sim 10^{4}$ K as predicted by
simulations [Sofue et al. 2016] (Fig. 3). This cut occurs at altitudes as low
as $|z|\sim<200$ pc, or $|b|<\sim 2^{\circ}$, which is difficult to detect
even it is present because of the strong absorption by the foreground disc gas
as above.

In either case, the observed apparent cut of soft X-ray brightness near the
galactic plane makes contrast to the continuous radio continuum ridge of the
NPS intersecting the disc.

The latitudinal variation of the X rays due to interstellar extinction has
been calculated using the hydrodynamical simulation as shown in Fig. 1. Fig. 2
compares the expected variation of intrinsic emissivity along the simulated
post shock (the densest inner shell of the simulation) by a thin line with
that after correction for the interstellar extinction by a thick line. Here,
the galactic disc is expressed by a plane parallel gas layer of thickness of
100 pc having absorption coefficient of $5\times 10^{-3}$ pc-1, or $\tau=1$ by
200 pc in the galactic plane. We plot the observed ROSAT intensities along the
NPS by grey dots, which are well fitted by the simulation with the disc
extinction.

Figure 2:  ROSAT X-ray count rates at 0.75 keV along the NPS in arbitrary unit
by dots (Snowden et al. 1997) compared with the intrinsic emissivity along the
post shock from simulation of the starburst model by the thin line (Sofue et
al. 2016), and with that after correction for the interstellar extinction due
to the galactic disc by the thick line. Compare the relative latitudinal
variations, as the absolute values are arbitrary.

From the strongly extinct nature near the Galactic plane as proved by Fig. 1
and 2, we may argue that the X-ray maps are not appropriate for studying the
physics of GSB interacting with the disc.

###  2.2 Temperature variation with height

On the other hand, X-ray maps at high latitudes may be used to represent and
study various physical properties of the GSB. Fig. 3 shows gaseous temperature
along the NPS plotted against the galactic latitude as observed with SUZAKU
[Kataoka et al. 2013, Tahara et al. 2015, Akita et al. 2018] and HaloSat
[LaRocca et al. 2020]. In all the data plotted here, the X-ray spectra are
modeled by a same three-component plasma model: APEC1 + Wabs $\times$ (APEC2 +
PL), where Wabs represents the Galactic absorption, APEC1 is an unabsorbed
thermal component that represents the Local Bubble (LB) emission, APEC2 is an
absorbed thermal component that represents the shocked halo gas like NPS, and
PL is the contribution from Cosmic X-ray background. Note that the temperature
of APEC2 is nearly constant at $kT\sim 0.3$ keV, or $T\sim 3\times 10^{6}$ K
(but see a slight different version of model fitting in which APEC2 is further
divided into $kT$ $\simeq 0.2$ and 0.4 keV components corresponding to
unshockd and shocked halo gass [Miller & Bregman 2016]). According to the
bipolar-hypershell model, the high temperature can be reasonably attributed to
heating of the halo gas by explosive event in the GC associated with a giant
shock wave expanding at $\sim 200-300$ km s-1[Sofue 1980, Sofue et al. 2016,
Sarkar 2019, Zhang & Guo 2020].

According to the simulation near the Galactic plane at $z<\sim 0.2$ kpc, the
temperature is expected to be significantly lower than that in the halo due to
the stronger deceleration of the shock wave inside the disc, causing decrease
in the temperature, and so the cooling rate and emissivity of X rays (figure
3). However, the current X-ray observations cannot measure the intrinsic
behavior at lower-latitudes than $|b|\leq 10^{\circ}$ because of the strong
absorption by the foreground galactic disc as in figure 2.

Figure 3: Gas temperature against latitude ($={\rm atan}(z/8{\rm kpc})$) at
the front shock (full line) and post shock (dashed line) from simulation of
the starburst model, compared with observed X-ray temperatures from SUZAKU
(circles) (Kataoka et al. 2013; Tahara et al. 2015; Akita et al. 2018) and
HaloSat (LaRocca et al. 2020).

As shown in Fig. 3, observed X-ray gas temperature of $kT$ $\simeq$ 0.3 keV is
surprisingly uniform in the halo against latitude, and is roughly consistent
with the ”post shock” temperature predicted from the simulation of the
starburst origin model.

On the other hand, temperature of ”front shock” is generally much higher. In
fact, MHD simulation predicts that the shock velocity gets even faster,
sometimes more than $\simeq$1000 km s-1, or $kT$ $\simeq$ 1 keV at high $b$
owing to substantial decrease of the halo gas density where the shock
propagates. Note that the cooling time scale of the X-ray emitting gas is
expressed as

$t_{\rm cool}\simeq 50\left(\frac{v_{sh}}{300~{}{\rm
km~{}s^{-1}}}\right)^{\alpha}\left(\frac{n}{0.01~{}\rm{cm^{-3}}}\right)^{-1}~{}\rm{Myr},$
(1)

for halo gas of subsolar metallicity $Z$ $\simeq$ 0.2 $Z_{\odot}$ [Draine
2011, Kataoka et al. 2021], where $\alpha=3.4$ at lower temperature than
$kT\sim 1$ keV or velocity $v_{\rm sh}\sim 1000$ km s-1 and $\alpha=1$ beyond.
Assuming $v_{sh}$ $\sim$ 1000 km s-1 and $n$ $\simeq$ 0.001 cm-3 for front
shock at high $b$, we obtain $t_{\rm cool}$ $\simeq$ 1.5 Gy, that is much
longer than dynamical scale for the bubbles’ formation. Hence the gas heated
by front shock can be effectively adiabatic during the expansion.

The emissivity of the X-ray emitting gas scales as $\propto$ $n^{2}$, thus
most of the observed emission comes from post-shock region rather than front
shock, where $kT$ $\simeq$ 0.3 keV plasma dominates. The cooling time is then
$\sim 12$ My, still long enough for the gas being regarded to be adiabatic.

The observed offset between radio and X-ray NPS is attributed to the position
of the front shock and post shock where the radio and X-ray emission becomes
strongest. In fact, our MHD simulation [Sofue et al. 2016] predicts $\sim$ 1
kpc offset at $z$= 4 kpc (or $b$ $\simeq$30∘), that corresponds to $\simeq$ 7
deg. This is roughly consistent with the observational radio/X-ray offset of
$\simeq$ 5 deg at the brightest part of the NPS [Kataoka et al. 2021].

###  2.3 Relieved X-ray maps and innermost X-ray spurs

In Fig. 5 we show a relieved map of soft X-ray intensity at 1 to 2.3 keV as
produced from the eROSITA all sky map [Predehl et al. 2020]. The relieving was
applied by shifting the image in the longitude direction by a couple of
angular resolution, then subtracting it from the original image. By this
simple procedure, the strong and relatively uniform galactic background is
subtracted and small scale structures are enhanced. Note, however, that the
images may not be used for quantitative measurement of the intensity, but is
useful only for morphological study of the edged structures.

Figure 4: Relieved X-ray image produced from eROSITA all sky map at 1 to 2.3
keV using the color image from (Predehl et al. 2020).

Figure 5: Horizontal-relieved map of $100^{\circ}\times 100^{\circ}$ region
around GC at 1420 MHz (Reich et al. 2001).

In the outer region at $|l|\sim>20^{\circ}$, X-ray spurs are coherently
aligned along the radio spurs with significant displacements toward inside the
bubbles [Kataoka et al. 2021]. Symmetric to the North Polar Spur (NPS) with
respect to the GC, a bright X-ray southern spur at $\sim$G$-30-20$, drawing a
clear arc, is associated with the South Polar Spur West (SPS-W) in radio
continuum. An X-ray spur at $\sim$G$-30+20$ is associated with the radio spur
NPS-W, and the X-ray spur at $\sim$G$+30-30$ corresponds to SPS-E in radio. As
a whole, both the radio and X-ray GSB compose an $\Omega$/  $\Omega$ shape,
making a symmetric dumbbell structure around the GC. Note, however, that they
do not draw an 8 shape crossing at the nucleus. However, low-latitude
extension below $|b|\sim<5^{\circ}$ of the X-ray spurs of the GSB and their
connection to the Galactic disc at $l\sim\pm 20-30^{\circ}$ cannot be studied
from the maps because of the strong absorption by the disc gas.

On the contrary, in a more central region at $|l|\sim<20^{\circ}$, we find
numerous vertical X-ray structures, which include the near-GC spurs at G+10-05
(SE Claw) and G-05+05 (NW Clump) [Kataoka et al. 2015] and the X-ray wind
horns [Bland-Hawthorn and Cohen 2003]. It is interesting to point out that
these vertical features are located inside the HI wind cone. High-velocity
ionized gas observed below the GC [Savage et al. 2017] could be related to
these near-GC X-ray spurs.

###  2.4 Origin of energy injection: AGN or starburst?

The GC explosion model of the large scale bubbles seems to be promising, but
origin of energy injection is yet to be understood, either by an AGN-like
outburst or from nuclear starburst in the GC. In general, AGN-like models
assumes that energy is released as a single point-like explosive event with a
short time-scale in the nucelus, whereas starburst assumes that the energy is
continuously released and injected to the GC over Myr timescale. While the
total amount of energy required is almost the same, AGN model predicts younger
bubble ages of 1$-$3 Myr owing to fast expansion velocity of 1000-10000 km s-1
[Guo & Mathews 2012]. Such high velocity is apparently inconsistent with the
observed X-ray temperature as discussed above, but situation can be
significantly moderated by assuming constant AGN activity or continuous energy
injection over Myr, leading to an expected plasma temperature of $kT$ $\simeq$
0.4 keV [Zhang & Guo 2020].

In this context, starburst models more naturally explain observed X-ray
temperature and overall morphology of bubbles owing to continuous star forming
activity over $\gtrsim$10 Myr [Sarkar et al. 2016, Sofue et al. 2016, Sarkar
2019]. In this scenario, the shock velocity and temperature in the highest
compression region in the post shock wave is much smaller than that at the
shock front by AGN models. The heated plasma can be gradually dissipated
and/or cooled to form quasi-isothermal bubbles as observed in X-rays. However,
it has been argued that the observed subsolar metallicity in the bubbles, $Z$
$\simeq$ 0.2 $Z_{\odot}$, seems to be difficult to be understood within the
framework of starburst models [Inoue et al. 2015].

It seems that difference between AGN and starburst models simply depends on
the energy injection time scale, rather than physical origin itself. Thus an
alternative idea could explain the isothermal nature by assuming recurrent and
continuous AGN activity as a source of energy injection over Myr time-scales.
In fact, such recurrent activity are often observed in number of radio
galaxies wherein multiple/mini radio lobes are found in quite different
spatial size (e.g., Centaurus A: [Israel 1998] for a review). The X-ray
chimney, which looks like mini-bubbles, recently discovered in the GC is also
collaborating this idea [Ponti et al. 2019]. It is however difficult to
clearly distinguish AGN and starburst scenario without detailed measurement of
physical conditions of X-ray emitting plasma, such as the gradient of
temperature, metallicity and ionization parameters against both latitude and
longitude directions.

##  3 Radio GSB intersecting the disc

###  3.1 Radio bubbles

Fig. 5 shows a relieved image of radio continuum emission at 1420 MHz for a
$100^{\circ}\times 100^{\circ}$ region centered on the GC made from the Bonn-
Stockert-Argentine all-sky survey [Reich et al. 2001]. Almost identical
features are observed in the relieved all-sky maps from the all-sky radio
surveys at 408 MHz [Haslam et al. 1982] and 2300 MHz [Jonas et al. 1998].

The north-eastern edge of the GSB, or the North Polar Spur (NPS), is visible
as the brightest radio and X-ray ridge near $l=30^{\circ},\ b=+20^{\circ}$
(hereafter, G30+20). The western counterpart to NPS (NPS-W) is recognized near
G-30+20 extending from the galactic plane to positive latitude. Their southern
counterparts are recognized around G30-20 and G-30-20, which were called the
South Polar Spurs East and West (SPS-E, SPS-W). These four spurs (NPS, NPS-W,
SPS-E and SPS-W) compose a bipolar double-horn structure symmetric around the
galactic plane and the rotation axis of the Galaxy.

###  3.2 Intersection with the galactic plane

The radio continuum spurs composing the GSB extend toward the galactic plane
without significant absorption, and merge with the complex disc composed of
numerous radio sources and extended structures. The northern radio spurs
apparently stop at the Galactic plane, and appear to be connected to the
southern spurs, although the details at $|b|\sim<1-2^{\circ}$ are not well
established in the figures mainly because of the contamination by the fore-
and background radio sources through the Galactic disc.

Among the spurs composing the GSB, the NPS is most prominent about the
relation to the Galactic disc. The extension of the radio NPS toward the
Galactic plane has been studied using background-filtered continuum maps
observed with the Effelsberg 100-m telescope [Sofue & Reich 1979, Sofue 2000].

In order to see the intersection of NPS with the galactic plane in more
detail, we show in Fig. 6 a relieved 2.8 GHz map made from the survey with the
100-m telescope at a resolution of $4^{\prime}.3$ [Reich et al. 1990] and
relieved 4.8 GHz map from the Urmqi 25-m telescope at a resolution of
$9^{\prime}.5$ [Sun et al. 2011]. In the bottom panel, we show linearly
polarized intensity map at 4.85 GHz [Sun et al. 2011], indicating that the
spur’s non-thermal nature of synchrotron origin. The polarization angle of the
E vector along the spur lie at PA$\sim 130^{\circ}$, indicating that the B
(magnetic) vector is at PA$\sim 40^{\circ}$, or the magnetic field is aligned
along the spur, if the Faraday rotation is negligible.

At very low latitudes of $|b|\leq\sim 1^{\circ}$, the radio continuum NPS,
both in intensity and polarization, shows no signature of extension across the
galactic plane. The fact that the NPS terminates at G21+01 at its brightest
and sharpest end indicates that the NPS has its origin in the Galactic disc.

Figure 6: [Top] 2.8 GHz relief map from Bonn-100 m telescope survey (Reich et
al. 1990) near the root of the NPS. The arrows points the extension of the
NPS. [Middle] 4.8 GHz relief from Urmqi 25-m telescope survey (Sun et al.
2010). [Bottom] Same, but polarized intensity. Arrows point the NPS ridge. A
bright source at $\sim$G19+03 is a foreground HII region, which is not visible
in polarization.

###  3.3 Magnetic fields

Besides the transparency for studying the morphological detail even through
the galactic disc, an advantage to use radio continuum emission is its direct
relation to the magnetic field in the shell.

If cosmic rays (CR) are produced by turbulent acceleration in the compression
site of magnetic fields, we may assume equipartition between CR and magnetic
energy densities. Then, the magnetic strength can be estimated using the
relation [Sofue et al. 2019],

$B\sim 3.1\times 10^{8}(\nu/1\ {\rm GHz})^{-1/7}\epsilon^{2/7}\ \mu{\rm G},$
(2)

where $\nu$ [in GHz] is the radio frequency and
$\epsilon=4\pi\int\Sigma_{\nu}d\nu/L\sim 4\pi\nu\Sigma_{\nu}/L$ [in erg cm-3
s-1] is the radio emissivity for line-of-sight depth $L$ of the emission
region. Here, $\Sigma_{\nu}=2k{T_{\rm B}}_{\nu}/\lambda^{2}$ is radio surface
brightness at $\nu$, ${T_{\rm B}}_{\nu}$ the brightness temperature, and
$\lambda$ the wavelength.

Inserting observed excess in brightness temperature of ${T_{\rm B}}\sim 0.5$ K
at 2.3 GHz [Jonas et al. 1998], $\sim 1.4$ K at 1.4 GHz [Reich et al. 2001],
and $\sim 40$ K at 408 MHz [Haslam et al. 1982] toward the NPS ridge at $b\sim
30^{\circ}$, we obtain volume emissivity of $\epsilon\sim 1.5\times 10^{-29}$,
$\sim 0.96\times 10^{-29}$, and $\sim 0.68\times 10^{-29}$ erg cm-3s-1,
respectively. These values yield magnetic strength of $B\sim 1.6$ $\mu$G for
an assumed line-of-sight depth of $L\sim 0.5$ kpc. Then, the magnetic energy
density, which is equal to CR energy density, is $u_{\rm mag}=u_{\rm CR}\sim
B^{2}/8\pi\sim 1.0\times 10^{-13}$ ergs cm-3. This is smaller than that of the
X-ray emitting thermal gas with $kT\sim 0.3$ keV and $n\sim 10^{-2}$ cm-3,
having $u_{\rm thermal}\sim 5\times 10^{-12}$ ergs cm-3. The decay time of CR
electrons emitting at 1.4 GHz is estimated to be $t_{\rm CR\ decay}\sim u_{\rm
CR}/\epsilon\sim 0.4$ Gy. This is longer than the dynamical time of the shock
wave, $t_{\rm dyn}\sim 10$ My, but shorter than the cooling time of the
thermal gas, $t_{\rm cool}\sim 1.5$ Gy.

It is noted that much more energy may be reserved in CR protons rather than
electrons, but the emission from CR protons are hard to be observed. In fact,
recent MHD simulation on the SNR shock suggests that energy density in the CR
protons is several orders of magnitude larger than that of CR electrons,
depending on the SNR age in Sedov phase [Urosević et al. 2018]. Moreover,
given that only a small fraction of cold protons are picked up into the CR
acceleration process, we should also consider hidden contribution of $cold$
protons as a reservoir of CR protons. Thus the substantial mismatch between CR
electrons, magnetic fields and X-ray emitting gas may suggest that most of
energy may be stored by protons in the front shock of the NPS and GSBs.

Morphological comparison between X-ray and radio spurs on the sky has posed an
interesting question about the displacement between the X-ray and radio edges,
such that the X-ray shell is located significantly inside the radio shell
[Kataoka et al. 2021]. The large displacement may imply either that the
frozen-in nature does not hold between the gas and magnetic fields in the GSB.
Alternatively, the acceleration of CR electrons is local and their life time
is comparable with the GSB’s evolutionary time scale, which is, however,
marginal in view of the life time estimation in section 4.3.

##  4 3-kpc crater as a footprint by GSB on the disc

###  4.1 GSB’s root enveloped by HI hole and 3-kpc crater

Figure 7:  X-ray image by eROSITA (Predehl et al. 2020) at 1 to 2.3 keV
compared with HI volume density ($n_{\rm HI}\propto T_{\rm B}$) at tangent
velocities, showing the cylindrical cavity coincident with the bubbles (Sofue
2017).

We first overview the spatial relation of the GSB to the HI hole and 3-kpc
crater in the Galactic disc and halo [Lockman 1984, Lockman & McClure-
Griffiths 2016, Sofue 2017]. In Fig. 7 we overlaid a map of the HI volume
density ($n_{\rm HI}\propto T_{\rm B}$) along the terminal velocity circle
made from GASS HI survey [McClure-Griffiths et al. 2009] on the eROSITA X-ray
map [Predehl et al. 2020]. The low-latitude X-ray spurs composing the GSB
approximately coincide in position with the inner wall of the HI hole (conical
cavity). The inner HI disk is as thin as $\sim 1^{\circ}$ inside 3 kpc, but
the disc suddenly becomes thicker and denser at $l\sim 22^{\circ}$, attaining
a thickness of $\sim 10^{\circ}$, exhibiting a crater structure. The HI wall
composes vertical spurs extending to $b\sim\pm 10^{\circ}$ (1 kpc), coinciding
also with the radio spurs. These facts agree with the idea that the NPS and
its counter spurs compose a dumbbell-shaped shock front, which swept away the
interior HI disk inside $\sim 3$ kpc.

###  4.2 Longitude-velocity diagrams

Fig. 8 shows HI line longitude-velocity (LV) diagrams made from the Parkes All
Sky Survey (GASS) of the HI line emission [McClure-Griffiths et al. 2009] at
$b=0^{\circ}$ and $0^{\circ}.6$. The 3-kpc expanding ring (HI crater’s wall)
is evident as the tilted ellipse of HI ridge in the LV diagram, which crosses
$l=0^{\circ}$ at positive and negative high velocities of $\pm 50-60$ km s-1.

Figure 8:  LV diagrams in the HI line at $b=0^{\circ}$ and $+0.6^{\circ}$ from
GASS survey (McClure-Griffiths et al. 2009), and CO line at $b=0^{\circ}$ and
$+0^{\circ}.5$ from Colombia survey (Dame et al. 2001) . The 3-kpc expanding
ring is indicated by the dashed ellipse. Figure 9: CO-line LV diagrams at
$b=+0^{\circ}$ (top) and $0^{\circ}.5$ (bottom) from the FUGIN survey (Umemoto
et al. 2017). Note the clear cut of terminal-velocity clouds at $l=22^{\circ}$
as indicated by the arrows.

The negative velocity (approaching) arm is recognized at $b\sim\pm
0^{\circ}.5-0^{\circ}.6$ both in CO and HI LVDs, while positive velocity arm
is fainter at $\pm 0^{\circ}.5$. In the off-plane regions at $b=\pm
0^{\circ}.5-0^{\circ}.6$ at $|l|\sim\leq\sim 20^{\circ}$, there appears a
vacant region of tangent-velocity gases, marked as HI hole. Such an LV
behavior implies that the gases are empty inside 3 kpc at these latitudes.

Fig. 9 shows LV diagrams of the ${}^{12}{\rm CO}(J=1-0)$emission at higher
resolution from the FUGIN222Four-receiver system Unbiased Galactic Imaging
survey with the Nobeyama 45-m telescope [Umemoto et al. 2017] CO survey around
the northern tangent point of the 3-kpc ring. The figure indicates that the
upper layer of the molecular disc is sharply truncated at $b\sim 0^{\circ}.5$
(65 pc).

###  4.3  Explosion vs bar origin

The non-circular motion exhibited by the elliptical feature in the LV diagram
in the 3-kpc arm has been interpreted in two ways.

One idea is that the ring is expanding, accelerated by an explosion at the GC
[Sanders & Prendergast 1974, Bania 1980, Sofue 1984]. The gas is approaching
to the Sun in the near side at negative radial velocity of $-53$ km s-1, while
receding in the far side at positive velocity at $+56$ km s-1[Dame & Thaddeus
2008]. This motion produces a simple LV ellipse, as shown by the dashed line
in figure 8, which reasonably traces the observed HI LV ellipse. The total
mass of the CO and HI rings is on the order of $\sim 10^{8}M_{\odot}$, and
their kinetic energy of the expanding motion at 50 km s-1is $\sim 5\times
10^{54}$ erg [Sofue 2017]. This requires only a small portion, $\sim 0.1$, of
the total energy to drive the ring’s expansion. The angular momentum problem
may be eased, if the ring is a focusing wave front and not an accumulated ring
[Sofue 1977].

Another idea is that the non-circular LV feature is produced by a bar
potential [Weinberg 1992] according to the various families of ”x” orbits
[Contopoulos 1956]. Simulations have shown complicated LV features, known as
the parallelogram, whose orientation depends on the position angle, mass,
axial lengths and ratio of the bar in a complex way [Binney et al. 1991,
Athanassoula & Bureau 1999, Sormani et al. 2015].

We may also combine the two ideas in such a way that the halo structures are
produced by the explosion, while the disc’s LV feature by the bar. The bar
promotes accretion of the disc gas to GC, and triggers starburst and/or AGN,
producing the GSB. However, in this case, we must assume 1) that the exact
coincidence of the 3 kpc arm with the roots of GSB and HI cylinder happened by
chance, 2) that the gigantic explosion did not give any imprint on the disc,
and 3) that the bar and flow orientation must be special to produce the
observed simple LV oval.

Considering these, particularly the fact that the observed 3-kpc LV oval is
fitted by a simple ellipse fairly well as in figure 8, not particularly
requiring a parallelogram, we here attempt to explain the GSB, HI cylinder,
3-kpc ring, and the HI+CO crater structure by a single hypothesis of the GC
explosion. Note that the explosion model explains the 3D structures in the
halo and disc, whereas the current bar simulations cannot touch upon the
vertical structures. A more sophisticated 3D model incorporating a bar and GC
blast would be desirable, although it is beyond the scope of this paper.

###  4.4 Density and disc thickness jumps at 3 kpc crater

Fig. 10 shows longitudinal variations of the HI and CO line brightness
temperatures $T_{\rm B}$ at the tangent-point velocities plotted against
longitude in the galactic plane and at $b=\pm 0^{\circ}.5-0^{\circ}.6$. We
recall that $T_{\rm B}$ is proportional to the mean volume density in the line
of sight depth through the tangent point gas as

$n_{i}=\frac{dN_{i}}{dx}=\frac{dN_{i}}{dv}\frac{dv}{dx}=X_{i}T_{\rm
B}\frac{dv}{dx},$ (3)

where $i$=HI or CO with $X_{i}$ being the conversion factors, $x$ and $v$ the
distance along the line of sight and radial velocity, respectively. The
velocity gradient is on the order of $dv/dx\sim 5$ km s-1kpc-1 for a velocity
gradient of $\sim 5$ km s-1and a line of sight depth across the tangent-
velocity gas of $\sim 1$ kpc. For HI with $X_{\rm HI}=1.8\times 10^{18}$ H
cm-3 [K km s-1]-1 and $T_{\rm B}\sim 80$ K, we obtain a mean density of
$n_{\rm HI}\sim 0.3$ H cm-3. For CO with $X_{\rm CO}\sim 2\times 10^{20}$ H2
cm-2 [K km s-1]-1 and $T_{\rm B}\sim 5$ K, we obtain $n_{\rm H_{2}}\sim 2$ H2
cm-3.

Figure 10: [Top] Longitudinal variation of $T_{\rm B}$ of HI gas (GASS) at
tangent points at different latitudes, $b=0^{\circ}$ and $\pm 0^{\circ}.6(\pm
78$ pc) . [Middle] FUGIN ${}^{12}{\rm CO}(J=1-0)$$T_{\rm B}$ at tangent points
at $b=+0^{\circ}.5$ (+65 pc). [Bottom] FUGIN ${}^{12}{\rm CO}(J=1-0)$intensity
at $V_{\rm lsr}\geq 100$ km s-1at $b=0^{\circ}$ and $0^{\circ}.5\ z=+65$ pc.
Note the jump of the gas density at $l\sim 22^{\circ}$ corresponding to the
3-kpc crater, within which the HI and CO gas discs are truncated.

The $T_{\rm B}$ plots in Fig. 10 show a step-like increase of the volume
densities of HI and CO gas from lower to higher longitudes at $l=22^{\circ}$.
In the Galactic plane ($b=0^{\circ}$), the density jump of the molecular gas
takes place from a finite but small value inside $22^{\circ}$ to higher
values. On the other hand, the off-plane molecular gas is almost empty inside
$=22^{\circ}$, but suddenly attains finite values beyond this longitude. Such
a density profile indicates that the molecular gas disc is thin inside 22∘,
and is surrounded by a thicker disc beyond the longitude, making a wall-like
distribution with a radius of 3 kpc.

Cross section of the gas distribution in the 3-kpc crater can be more directly
displayed by channel maps in the HI and CO-line emissions at the terminal
velocity. Fig. 11 shows $I_{\rm CO}$ maps integrated at $V_{\rm lsr}\geq 60$
km s-1from the Columbia CO survey, and the same but at $V_{\rm lsr}\geq 100$
km s-1from the Nobeyama FUGIN survey around $l\sim 22^{\circ}$. A sudden jump
of the thickness of CO disc is recognized at $l=22^{\circ}$ such that the
thickness by FWHM (full width at half maximum) of $\sim 0^{\circ}.5$ at $l\leq
22^{\circ}$ to $\sim 1^{\circ}.5$ at $l>22^{\circ}$.

The sudden change of the thickness revealed in the FUGIN CO map suggests that
the upper layer of the molecular disc has been sharply truncated by sweeping
mechanism inside $l\leq 22^{\circ}$ due to a shock wave from the GC or a
horizontal wind blowing in the halo.

Figure 11:  [Top] ${}^{12}{\rm CO}(J=1-0)$$I_{\rm CO}$ maps integrated at
$V_{\rm lsr}\geq 60$ km s-1from Columbia survey (top), and [Bottom] $I_{\rm
CO}$ from FUGIN survey integrated at $V_{\rm lsr}\geq 100$ km s-1. Contours
are at 10, 30 and 50 K km s-1. Note the step-like increase of the disc
thickness at $l\sim 22^{\circ}$ as traced by the lowest contours, indicating
truncation of the upper layer inside 3 kpc.

The thickness variation may be more clearly seen by comparing vertical cross
sections at different longitudes outside and inside of the 3-kpc arm at G22.
Fig. 12 shows the cross section at G22.5 compared with that at G21. The former
shows a wide extent of the disk, having an outskirt to high latitudes. On the
other hand, the cross section at G21 shows a plateau-like distribution,
sharply cut at $\pm 0^{\circ}.25$ from the centroid, suggesting truncation of
the outskirt.

Figure 12: Latitude cross sections of $I_{\rm CO}$ ($V_{\rm lsr}\geq 60$ km
s-1) at G21 and G22.5 around the 3-kpc ring from FUGIN, showing drastic change
of the vertical extent of the molecular disc from inside to outside the 3 kpc
ring at $l=22^{\circ}$.

Thus, the horizontal (longitudinal) and vertical (latitudinal) variations show
that the density and thickness of the galactic gas disk, both in HI and CO,
increase suddenly at $l\sim 22^{\circ}$. The variation is more step-like in
the off-plane region at $|b|\sim 0^{\circ}.5$ than in the galactic plane.

From Fig. 12 the molecular disk thickness inside G22 is as thin as $\pm
0^{\circ}.25$ (30 pc) and increases to $\sim 1^{\circ}$ outside $22^{\circ}$.
This indicates that the galactic disc of neutral gas (HI and CO) has a crater
structure, having a hole inside G22 above the disk at $|z|>\sim 30$ pc
($0^{\circ}.25$). As shown in Fig. 10, the HI hole’s wall has a radius of
$r=3$ kpc and thickness $\sim 260$ pc (rising full width of $\sim 2^{\circ}$
measured at $b=0^{\circ}.6$). On the other hand, the CO hole has a much
sharper wall with thickness of $\sim 30$ pc ($\sim 0^{\circ}.2$ at
$b=0^{\circ}.5$).

###  4.5 Formation of the crater and energetics

It has been shown that the root of the GSB can be traced by radio continuum
ridges intersecting the Galactic plane. The NPS root was shown to exactly
coincide with the 3-kpc expanding ring composed of HI and molecular gases.
Fig. 13 illustrates a schematic view of the interaction of the GSB with the HI
and CO disc based on the observed data shown in the previous subsections. In
the upper halo, the GSB makes giant shells with the radio continuum being
located outside X-rays [Kataoka et al. 2021], and is enveloping the
$\gamma$-ray bubble [Su et al. 2010]. The upper layer of the galactic disc is
swept away and accumulated in the wall of the 3-kpc crater, creating the HI
and CO rings. By the focusing effect of the shock wave on the disc [Sofue
1977, Sofue 2020], the front near the disc makes a compressed ring with the
expansion velocity at $\sim 50$ km s-1. However, the thin and densest area
near the Galactic plane at $|b|\lesssim 1^{\circ}$ remains almost unperturbed
because of the rapid deceleration of the shock wave in the plane inside 1-2
kpc from the GC.

Figure 13: Schematic illustration of the relation of GSB (NPS) to the GC and
3-kpc expanding ring (Norma Arm) and crater.

Using the observed brightness temperatures and three-dimensional distribution
expected from the apparent extent on the sky and assumed distance, the mass of
the 3-kpc ring has been estimated to be $1.5\times 10^{8}M_{\odot}$ [Sofue
2017]. The corresponding kinetic energy is estimated to be on the order of
$\sim 5\times 10^{54}$ ergs for the expanding motion at 50 km s-1[Sofue 2017].
This energy is an order of magnitude smaller than the currently estimated
total energy required for the formation of the GSB of $\sim 10^{55-56}$ ergs.
Hence, energetically, it is reasonable to assume a common origin for the
formation of the GSB and the 3-kpc expanding ring with the HI and CO crater in
the galactic disc.

##  5 Discussion

###  5.1 Unified view of the GSB in X, radio, HI and CO

We have re-visited the GSB using the all-sky X-ray and radio continuum surveys
as well as HI and the newest CO-line data of the galactic disc. From the wide
coverage of the observed energy ranges as well as by the spectral line
information, we are now able to obtain a unified view of the GSBs as follows.

The X-ray data are crucial to derive thermal and kinematic properties of the
GSB as a shock wave from the GC propagating through the halo gas. Gaseous
temperature can be used to calculate the shock wave velocity $v$, which is
related to the injection energy $E_{0}$ at the GC and accumulated total mass
$M$ of the shell by $E_{0}\sim 1/2Mv^{2}$ for an adiabatic shock wave as
reasonably assumed from the sufficiently long cooling time. From X-ray
spectroscopy, we know that $kT\sim 0.3$ keV (Fig. 3), and the intensity yields
thermal electron density on the order of $n\sim 10^{-2}$ cm-3 for assumed
radius $\sim 4$ kpc and thickness $\sim 1$ kpc. These estimations lead to
$M\sim 2\times 10^{8}M_{\odot}$ and $E_{0}\sim 4\times 10^{55}$ ergs. Under
the adiabatic assumption, the age of the GSB is approximately estimated from
the current expansion velocity as $t\sim 0.6r/v\sim 11$ My.

Although X-ray data are useful for discussing the shock wave propagation and
energetics, they strongly suffer from the interstellar extinction by heavy
elements in the galactic disc, where the optical depth steeply increases
toward the galactic plane as $\propto{\rm sec}\ |b|$. This results in the dark
absorption belt along the galactic plane in the all-sky X-ray maps, making it
difficult to obtain information about the interaction of the high-temperature
gas with the galactic disc.

Such difficulty can be eased by observing radio continuum emission that can be
safely assumed to be transparent in the whole Galaxy including the disc. We
have thus showed that the roots of the GSB are connected to the tangential
directions of the 3-kpc expanding ring as revealed in the HI and CO line
emissions.

The radio ridge of the NPS was shown to be extending straightly toward the
galactic plane, getting sharper and brighter and terminating at $b\sim
1^{\circ}$ after merging into the bright disc emission. From the intersection
of the GSBs’ spurs at longitudes of $\sim\pm 20^{\circ}$, we showed that the
whole GSB structure composes an $\Omega$/  $\Omega$ shape.

The distribution of HI and molecular gases along the tangent-velocity circle
showed that the ring composes a crater-like structure in the disc around the
GC. The 3-kpc ring and crater wall are expanding at $\sim 50$ km s-1as
inferred from the LV diagrams. Kinetic energy of the expanding motion of the
ring is estimated to be on the order of $5\times 10^{54}$ ergs. Namely, a
small portion of the total energy required to drive the whole GSB
($10^{55-56}$ ergs) is sufficient to drive the expanding ring and to create
the crater.

We emphasize that the GC explosion model explains the coherent 3D structures
in the disc and halo around 3 kpc, which include the sudden increase in the
disc thickness at $l\sim 22^{\circ}$ beyond horizontal truncation of the upper
layer inside 3 kpc and the cylindrical HI wall around the central cavity as
well as the non-circular kinematics revealed in the LV digrams.

###  5.2 Galactic-scale feedback

Finally, we will comment on the implications of the interaction between the
GSB and the galactic disc, assuming that the expanding ring of the gas disk is
due to a GC explosion. Current research of GSB in the literature has been
devoted to discussing its impact on the galactic halo. However, in the
following we highlight the impact on the disc and its activity. This may
somehow akin to discussing the impact of an explosion on the earth not only on
the atmosphere, but also on structures on the ground.

The shock wave velocity is variable with the gaseous density in such a way
that the velocity decreases with the density. Accordingly, the ray path of a
wave is refracted due to varying propagation velocity in the disc and halo.
This causes a converging flow of the shock front toward the galactic disc,
making a ring of focusing waves [Sofue 2020].

If the wave amplitude is small enough, or it is a linear wave like sound and
MHD waves, the focal length $f$ is approximately related to the scale height
of the disc $h$ as $f\sim 2.2\ h$ [Sofue 1977, Sofue 2020]. Thus, besides the
large-scale energetic flow as huge bubbles into the halo, some fraction of the
released energy at the GC propagates the galactic disc and lower halo,
repeating focusing onto rings of radii $r_{i}\sim 2.2\ i\ h$ with $i=1,2,...$.

The 3-kpc ring could be one of such focal rings driven by the GC activity,
although its amplitude is non-linear and the simple focusing law may not apply
and the focal length would be longer than $2.2\ h$. Also, after passing a
focal ring, the wave diverges again into the disc and halo, propagates further
outward, repeating focusing.

It is, therefore, expected that some fraction of the released energy at the GC
by starburst or AGN activity is efficiently confined in the disc and lower
halo, and repeat convergence onto disc-wide focal rings. We emphasize that,
besides the GSB’s effects on the halo and intergalactic space, the galactic-
scale feedback of the GC activities back onto the disc is also important for
considering the structure and evolution of the Galaxy [Sofue 2020]. For
example, convergence of the released energy from the GC, even though a
fraction of the total amounting to $\sim 10^{55-56}$ ergs per $\sim 10$ My,
would cause various dynamical effects in the disc, which include triggering of
star formation by compression in focal rings, and generation of the
interstellar turbulence.

##  6 Summary

We revisited the GC super bubbles (GSBs), being inspired by the eROSITA all-
sky X-ray observations [Predehl et al. 2020], and discussed their origin as
due to an explosion in the Galactic Centre. With a particular emphasis of the
interaction of the bubbles with the Galactic disc, we investigated the radio
continuum, X rays, HI surveys and the newest high-resolution CO-line survey of
the Galactic disc. We showed that the roots of the bubbles coincide with the
3-kpc expanding ring and crater in the HI and CO discs of the Galaxy. The thus
revealed 3D structure of the crater is explained by sweeping of the disc gas
by the GSB as a footprint of the explosive event in the Galactic Centre. The
origin of the GSB can be well modeled by a $\Omega$/  $\Omega$ shaped shock
wave from the GC with an explosive energy input of several $10^{55}$ ergs
$\sim 10$ My ago.

## Acknowledgements

The authors are indebted to Prof. A. Habe for the hydrodynamical simulation
for Fig. 1 as reproduced from our paper in 2016. J.K. acknowledges the support
from JSPS KAKENHI Grant Number JP20K20923.

## Data availabivity

The CO line data were taken from the FUGIN CO-line survey at url: https://nro-
fugin.github.io/. The eROSITA X-ray image was taken from the url:
https://www.mpe.mpg.de/7461761/news20200619.

## References

  * [Akita et al. 2018] Akita. M., Kataoka, J., Arimoto, M. et al. 2018, ApJ, 862, 88. doi:10.3847/1538-4357/aacd08
  * [Athanassoula & Bureau 1999] Athanassoula, E. & Bureau, M. 1999, ApJ, 522, 699. doi:10.1086/307677
  * [Bania 1980] Bania, T. M. 1980, ApJ, 242, 95. doi:10.1086/158446
  * [Binney et al. 1991] Binney J., Gerhard O. E., Stark A. A., Bally J., Uchida K. I., 1991, MNRAS, 252, 210
  * [Bland-Hawthorn and Cohen 2003] Bland-Hawthorn, J., and Cohen, M. 2003, ApJ, 582, 246
  * [Carretti et al. 2013] Carretti E., Crocker R. M., Staveley-Smith L. et al., Nature, 493, 66
  * [Cohen & Davies 1976] Cohen R. J., Davies R. D., 1976, MNRAS, 175, 1
  * [Contopoulos 1956] Contopoulos, G. 1956, ApJ, 124, 643. doi:10.1086/146271
  * [Crocker et al. 2015] Crocker R. M., Bicknell G. V., Taylor A. M., Carretti E., 2015, ApJ, 808, 107
  * [Dame et al. 2001 ] Dame, T. M., Hartman, D., Thaddeus, P. 2001, ApJ 547, 792.
  * [Dame & Thaddeus 2008 ] Dame, T. M. & Thaddeus, P. 2008, ApJL, 683, L143. doi:10.1086/591669
  * [Draine 2011] Draine, B. T., Physics of the interstellar and intergalactic medium, 2011
  * [García et al. 2014] García P., Bronfman L., Nyman L., Dame T. M., Luna A., 2014, ApJS, 212, 2
  * [Guo & Mathews 2012] Guo, F., & Mathews, W. G., 2012, ApJ, 756, 181, doi: 10.1088/0004-637X/756/2/181
  * [Haslam et al. 1982] Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, AAS, 47, 1
  * [Heywood et al.2019] Heywood I., et al., 2019, Natur, 573, 235
  * [Inoue et al. 2015] Inoue, Y., Nakashima, S., Tahara, M., et al. 2015, PASJ, 67, 56, doi: 10.1093/pasj/psv032
  * [Israel 1998] Israel. F. P., The Astrnomy and Astrophysics Review, 1998, 8. 237, doi: 10.1007/s001590050011
  * [Jonas et al. 1998] Jonas, J. L., Baart, E. E., & Nicolson, G. D. 1998, MNRAS, 297, 977. doi:10.1046/j.1365-8711.1998.01367.x
  * [Kalberla et al. 2005] Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, AA, 440, 775
  * [Kataoka et al. 2013] Kataoka, J., Tahara, M., Totani, T., et al. 2013, ApJ, 779, 57, doi:10.1088/0004-637X/779/1/57
  * [Kataoka et al. 2015] Kataoka, J., Tahara, M., Totani, T., et al. 2015, ApJ, 807, 77. doi:10.1088/0004-637X/807/1/77
  * [Kataoka et al. 2018] Kataoka J., Sofue Y., Inoue Y., Akita M., Nakashima S., Totani T., 2018, Galax, 6, 27
  * [Kataoka et al. 2021] Kataoka J., Yamamoto, M., Nakamura, Y., et al., 2021, ApJ, vol.908, 14, doi:10.3847/1538-4357/abdb31
  * [Lallement et al. 2016] Lallement, R., Snowden, S., Kuntz, K. D., et al. 2016, AA, 595, A131. doi:10.1051/0004-6361/201629453
  * [LaRocca et al. 2020] LaRocca, D. M., Kaaret, P., Kuntz, K. D., et al. 2020, ApJ, 904, L54. doi:10.3847/1538-4357/abbdfd
  * [Lockman 1984] Lockman F. J., 1984, ApJ, 283, 90
  * [Lockman & McClure-Griffiths 2016] Lockman F. J., McClure-Griffiths N. M., 2016, ApJ, 826, 215
  * [McClure-Griffiths et al. 2009] McClure-Griffiths N. M., et al., 2009, ApJS, 181, 398
  * [Miller & Bregman 2016] Miller, J. M. & Bregman, J. N., 2016, ApJ, 829, 9
  * [Nakashima et al. 2013] Nakashima, S., Nobukawa, M., Uchida, H. et al. 2013, ApJ, 773, 20
  * [Planck collaboration 2013] Planck collaboration, 2013, A$\&$A, 554, 139
  * [Ponti et al. 2019] Ponti, G., Hofmann, F., Churazoz, E. et al. 2019, Nature, 567, 347
  * [Ponti et al. 2021] Ponti, G., Morris, M. R., Churazoz, E. et al. 2021, A$\&$A, in press (arXiv:2101.05284)
  * [Predehl et al. 2020] Predehl, P., Sunyaev, R. A., Becker, W., et al. 2020, Nat., 588, 227. doi:10.1038/s41586-020-2979-0
  * [Reich et al. 1990] Reich, W., Fuerst, E., Reich, P., & Reif, K. 1990, AA Suppl. 85, 633
  * [Reich et al. 2001] Reich, P., Testori, J. C., & Reich, W. 2001, AA, 376, 861. doi:10.1051/0004-6361:20011000
  * [Sanders & Prendergast 1974] Sanders R. H., Prendergast K. H., 1974, ApJ, 188, 489
  * [Sarkar et al. 2016] Sarkar K. C., Nath B. B., Sharma P., Shchekinov Y., 2016, ApJ, 818, L24
  * [Sarkar 2019] Sarkar, K. C. 2019, MNRAS, 482, 4813
  * [Savage et al. 2017] Savage, B. D., Kim, T.-S., Fox, A. J., et al. 2017, ApJS, 232, 25. doi:10.3847/1538-4365/aa8f4c
  * [Snowden et al. 1997] Snowden S. L., et al., 1997, ApJ, 485, 125
  * [Sofue 1977] Sofue, Y. 1977, AA 60, 327.
  * [Sofue 1980] Sofue, Y. 1980, PASJ, 32, 79
  * [Sofue 1984] Sofue Y., 1984, PASJ, 36, 539
  * [Sofue & Handa 1984] Sofue, Y. & Handa, T. 1984, Nat., 310, 568. doi:10.1038/310568a0
  * [Sofue 1994] Sofue, Y. 1994, ApJ.L., 431, L91
  * [Sofue 2000] Sofue, Y. 2000, ApJ, 540, 224
  * [Sofue 2015] Sofue Y., 2015, MNRAS, 447, 3824
  * [Sofue 2017] Sofue, Y. 2017, PASJ, 69, L8. doi:10.1093/pasj/psx067
  * [Sofue et al. 2016] Sofue Y., Habe A., Kataoka J., Totani T., Inoue Y., Nakashima S., Matsui H., Akita M., 2016, MNRAS, 459, 108
  * [Sofue et al. 2019] Sofue, Y., Nakanishi, H., & Ichiki, K. 2019, MNRAS, 485, 924. doi:10.1093/mnras/stz407
  * [Sofue 2020] Sofue, Y. 2020, MNRAS, 498, 1335. doi:10.1093/mnras/staa2389
  * [Sofue & Reich 1979] Sofue, Y. and Reich, W. 1979 AAS 38, 251
  * [Sormani et al. 2015] Sormani, M. C., Binney, J., & Magorrian, J. 2015, MNRAS, 454, 1818. doi:10.1093/mnras/stv2067
  * [Su et al. 2010] Su, M., Slatyer, T. R., and Finkbeiner, D. P. 2010, ApJ, 724, 1044
  * [Sun et al. 2011] Sun X. H., Reich W., Han J. L., Reich P., Wielebinski R., Wang C., Müller P., 2011, A&A, 527, A74
  * [Tahara et al. 2015] Tahara, M., Kataoka, J., Takeuchi, Y. et al. 2015, ApJ, 802, 91, doi: 10.1088/0004-637X/802/2/91
  * [Umemoto et al. 2017] Umemoto, T., et al., 2017, PASJ, 69, 78
  * [Urosević et al. 2018] Urosević,D., Pavlović, M. Z. & Arbutina, B., ApJ, 2018, 855, 59, doi:10.3847/1538-4357/aaac2d
  * [Weinberg 1992] Weinberg, M. D. 1992, ApJ, 384, 81. doi:10.1086/170853
  * [Zhang & Guo 2020] Zhang, R. & Guo. F., ApJ, 2020, 894, 117, doi: 10.3847/1538-4357/ab8bd0