# Implementation of nonlocal non-Fourier heat transfer for semiconductor
nanostructures

Roya Baratifarimani  Zahra Shomali  Department of Physics, Basic Sciences
Faculty, Tarbiat Modares University, Tehran, Iran<EMAIL_ADDRESS>
###### Abstract

The study of heat transport in micro/nanoscale structures due to their
application, especially in Nanoelectronics, is a matter of interest. In other
words, the precise simulation of the temperature distribution inside the
transistors is consequential in designing and building more reliable devices
reaching lower maximum temperatures during the operation. The present study
constitutes a framework for micro/nanoscale heat transport study which leads
to the calculation of accurate temperature/heat flux profiles with low
computational cost. The newly non-dimensional parameter $\gamma$, presenting
the strength of the nonlocality, is utilized through the nonlocal DPL modeling
(NDPL). Alongside the calculating nonlocality coefficient, the factors also
appearing in DPL, including the temperature jump, phase lagging ratio, are
revisited. The factor $\gamma$ is found to have a linear relationship with
Knudsen (Kn) number, being 3.5 for Kn=10 and 0.035 for Kn=0.1. Although the
nonlocality is bold for the large Knudsen numbers, it also plays a vital role
for low Knudsen number structures especially at earlier times. Further, It is
obtained that intruding $\gamma$ is critical for obtaining accurate
temperature and heat flux distributions which are very close to the practical
results of Phonon Boltzmann equation.

###### keywords:

Dual Phase Lag , nonlocality , Phonon Boltzmann equation , Nanoscale heat
transport , Thermal management

##  1 Introduction

Micro/nanoscale heat conduction in nano-materials, especially the
semiconductors, has attracted significant attention in recent years. This
concern on one hand, relates to the practical application in nano-electric [1,
2], and on the other hand, provides worthwhile understanding of the
fundamental principles behind the nanoscopic heat carries’ behavior [3, 4].
The Fourier’s law, $\vec{q}(\vec{r},t)=-k\vec{\nabla}T(\vec{r},t)$, which is
established as the accurate model to simulate heat conduction in macroscopic
structures, fails when length/time scales are, respectively, comparable to the
phonon mean free paths and relaxation times [5, 6, 7].

Previously, many studies have been performed focusing on the heat transfer
temporal behavior at the nanoscale. Specifically, the constitutive macroscopic
equation introduced as the dual-phase-lag (DPL) has been suggested:

$\vec{q}(\vec{r},t+\tau_{\\_}q)=-k\vec{\nabla}T(\vec{r},t+\tau_{\\_}t)$ (1)

In Eq. 1, two phase lags of $\tau_{\\_}T$ and $\tau_{\\_}q$ represent,
respectively, the phase lag of the heat flux and the temperature gradient. In
comparison to the existing microscopic models, the lagging behavior is used in
the place of electron/phonon coupling in metals [8], umklapp and normal phonon
scattering [9], and other relaxation of internal energy [10]. Both thermal
phase-lag values are determined experimentally [11] and theoretically [12].

Meanwhile, significantly fewer researches were devoted to the spatial attitude
of the non-Fourier heat transfer. In principle, the quasi-ballistic heat flux
in a given location and at a given time intrinsically depends on the
temperature gradient in other places, in addition to the earlier times [13].
In other words, the constitutive law is no longer localized in non-diffusive
transport regime, but instead possesses “memory” in space alongside the time
memory. From a microscopic point of view, a convolution kernel $\kappa^{*}$,
which is called the nonlocal thermal conductivity, fully contains the spatio-
temporal memory of the heat flux concerning the temperature gradient [14].
This parameter appears in the following postulated constitutive law, or the
so-called flux-gradient relation (FGR), which holds when the conventional
Fourier law no longer does hold:

$q(x,t)=-\int_{\\_}{0}^{t}dt^{\prime}\int_{\\_}{-\infty}^{\infty}dx^{\prime}\kappa^{*}(x-x^{\prime},t-t^{\prime})\frac{\partial
T}{\partial x}(x^{\prime},t^{\prime})$ (2)

Here, the convolution kernel $\kappa^{*}$, is the “nonlocal thermal
conductivity kernel” of the medium. $\kappa^{*}$ proceeds as an intrinsic
material characteristic of the thermal conducting medium. Under the relaxation
time approximation, the Boltzmann transport equation conventionally follows
the postulated constitutive law. After Fourier and Laplace transformations and
inserting the analytical single pulse response of the RTA-BTE, one yields a
generic expression regarding the underlying microscopic phonon dynamics [14].
Once the kernel is obtained, the nonlocal heat transport is maintained. For
the thermal transport through the length and time scales that are relatively
larger than the phonon mean free paths and relaxation times, $\kappa^{*}$
reaches the constant value and, consequently, the FGR depicted by Eq. 2,
becomes localized and reduces to Fourier’s law.

On the other side, from a macroscopic viewpoint, the equivalency between the
spatial nonlocal behavior and the lagging response in time, had been confirmed
[15, 16]. However, due to the absence of a constitutive model supporting, both
were not considered concurrently. In 2010, Tzou and Guo, extended the concept
of thermal lagging, in time, to the nonlocal response in space [17]. The study
was endorsed under the phonon scattering, the Guyer–Krumhansl model, and the
thermomass model, which expresses the same nonlocal behavior as the first- and
second-order effects. At first, nonlocal behavior with thermal lagging was
motivated by the famous Guyer–Krumhansl model [9]. The model described via the
following equation, includes the parameters which characterize the heat wave
transition through the dielectric crystals:

$\textbf{q}+\tau_{\\_}R\frac{\partial\textbf{q}}{\partial
t}+(\frac{c^{2}\tau_{\\_}Rc^{2}}{3})\nabla
T=\frac{\tau_{\\_}R\tau_{\\_}Nc^{2}}{5}[\nabla^{2}\textbf{q}+2\nabla(\nabla.\textbf{q})].$
(3)

The above constitutive equation for heat transport in phonon systems,
introduces two relaxation times, $\tau_{\\_}R$, and $\tau_{\\_}N$, which,
respectively, describe the umklapp (momentum-nonconserving), and normal
(momentum-nonconserving) processes in a phonon framework. The relaxation times
are related to two phase lags of DPL via $\tau_{\\_}q$=$\tau_{\\_}R$ and
$\tau_{\\_}T$=9$\tau_{\\_}N$/5 [16]. The $\tau_{\\_}R\tau_{\\_}Nc^{2}$/5 on
the right-side of the Eq. 3, has a length squared dimension. The lengths
$\tau_{\\_}Rc$ and $\tau_{\\_}Nc$ are proportional to the phonon path traveled
during $\tau_{\\_}R$ and $\tau_{\\_}N$, and so
$\tau_{\\_}R\tau_{\\_}Nc^{2}$/5$\approx l^{2}$. The first-order effect in
”$l$” does not appear in Guyer–Krumhansl equation. It is necessary to conform
this spatial accomplishment together with the temporal effect to clear out the
heat transport in phonon systems.

Besides, the thermomass (TM) model [18] is also another framework that
supports the presence of the nonlocal behavior together with thermal lagging
while investigating heat transport. Through the TM model, where a finite mass
calculated from Einstein’s mass-energy relation is attributed to phonons, the
energy and constitutive equation in heat transport are derived, subsequently,
from the continuity and the momentum equation. In consequence, the continuity
equation reduces to the energy equation in the phonon gas:

$\frac{\partial\rho_{\\_}h}{\partial t}+\frac{\partial}{\partial x}(\rho
u_{\\_}h)=0\ \implies-\frac{\partial q}{\partial x}=C\frac{\partial
T}{\partial t}.$ (4)

Also, the constitutive equation in one-dimensional reads out:

$\rho_{\\_}{h}\overbrace{(\frac{\partial u_{\\_}h}{\partial
t}+u_{\\_}h\frac{\partial u_{\\_}h}{\partial
x})}^{\frac{Du_{\\_}h}{Dt}}+\frac{\partial P_{\\_}h}{\partial
x}+f_{\\_}h=0\;\implies\;\tau_{\\_}{TM}\frac{\partial q}{\partial t}-\\\
\Lambda C\frac{\partial T}{\partial t}+\Lambda\frac{\partial q}{\partial
x}-M^{2}k\frac{\partial T}{\partial x}+k\frac{\partial T}{\partial x}+q=0.$
(5)

Here, $\tau_{\\_}{TM}$, the lagging time in the thermomass model, is about two
orders of magnitude larger than the heat flux phase lag in the CV model. Also,
$l=qk_{\\_}p/[2\gamma C(CT)^{2}]$ and M are, respectively, the length
parameter, and the drift velocity thermal Mach number with respect to the
thermal wave speed. When the heat flux is replaced from Eq. 4 to Eq. 5, one
obtains an energy equation that contains only the temperature:

$\frac{k(1-M^{2})}{C\tau_{\\_}{TM}}\frac{\partial^{2}T}{\partial
x^{2}}=\frac{1}{\tau_{\\_}{TM}}\frac{\partial T}{\partial
t}+\frac{2l}{\tau_{\\_}{TM}}\frac{\partial^{2}T}{\partial x\partial
t}+\frac{\partial^{2}T}{\partial t^{2}}.$ (6)

The parameters involved in the above equations are $\rho_{\\_}h$, $u_{\\_}h$,
$p_{\\_}h$, and $f_{\\_}h$, which are, respectively, defined as the phonon gas
thermomass density, the drift velocity, the pressure of the phonon gas, and
the resistance force per unit volume. Although the wave behavior is understood
by the second order time derivative, $\partial^{2}T/\partial t^{2}$, Eq. 6
also consists of a mixed-derivative term, $\partial^{2}T/\partial x\partial
t$. This part is new to the field of microscale heat transfer. The standard
DPL model, considering only the phase lags, does contain similar mixed-
derivative terms, which are even orders with respect to $x$. As it is seen,
the length parameter, $l$, in Eq. 6, is proportional to the heat flux, $q$.

The presence of the length parameter in Eqs. 3 and 6 is explanatory for
including the nonlocal attitude in space, as well as the lagging time in the
CV model. Such nonlocality and lagging are implied in Fourier’s law through
the succeeding statement:

$q(x+\lambda_{\\_}q,t+\tau_{\\_}q)=-K\frac{\partial T}{\partial x}(x,t).$ (7)

The equation supports that the temperature gradient through the material
volume on position r and at time t, is proportional to the heat flux vector
flowing across the other volume element positioned at ($r+L$) at a later time
(t + $\tau_{\\_}q$). The parameters, $K$=$k(1-M^{2})$,
$\tau_{\\_}q$=$\tau_{\\_}{TM}$, and $\lambda_{\\_}q$=2$l$, ascertain that the
nonlocal model is equivalent to the TM model. In the next step, the
nonlocality is implied to more accurate DPL model [19, 20]. Consequently, the
Nonlocal DPL equation becomes:

$q(x+\lambda_{\\_}q,t+\tau_{\\_}q)=-K\frac{\partial T}{\partial
x}(x,t+\tau_{\\_}T).$ (8)

On the other hand, the heat transport study of the transistors is particularly
significant as their reliability is determined via the obtained maximum
temperature during the operation [21, 22, 23]. When the temperature and heat
flux distributions are calculated, we get information about how much
temperature increase, each position inside the single transistor or on a whole
die, experiences [24, 25, 26, 27, 28, 29, 30]. This temperature augmentation
in a MOSFET occurs due to the response to the existent self-heating [31, 32,
33, 34]. Consequently, predicting the temperature and heat flux distribution
preciously forasmuch as the reliability is determined can help the engineers
to design MOSFETs with optimal thermal conditions. Here, the nonlocal DPL
model has been utilized for the thermal investigation of the 1-D MOSFETs.
Newly nonlocality, phase lagging and, temperature jump coefficients are
derived. Explicitly, the obtained scaled parameters are found to be almost
different from that of the localized standard DPL model. The following study
confirms that the implication of the nonlocal effects notably makes the
simulation yield a more accurate results [35]. To be more precise, although
for low Knudsen number devices, the phenomenological DPL model brings out less
meticulous results, taking nonlocality in space for DPL, leads to much more
accurate temperature profiles and heat flux plots. In this paper, first, in
Sec. 2, the studied geometry and the relevant boundary conditions are given.
Then in Sec. 3, the mathematical modeling for NDPL is developed. Section. 4 is
devoted to explaining about the numerical method. Finally, the results are
pointed out in Sec. 5 and, the paper is concluded in Sec. 6.

##  2 Geometry and Boundary conditions

A one-dimensional MOSFET device is modeled in two cases. First, as shown in
Fig. 1 (a), a thin slab of silicon transistor with a constant initial
temperature of $T_{\\_}0=300\ K$ is studied. The system is the same as that of
the one investigated by Ghazanfarian and Abbassi [35]. The temperature of the
bottom boundary has abruptly increased. This results in the rearrangement of
the temperature distribution in the slab. The heat transport in the thin film
is studied considering the lagging responses and the nonlocality in space
simultaneously. During the calculation, the applied temperature at the left
boundary is kept constant at room temperature. The size of the transistor, L,
is selected as the characteristic length of the simulation and manages the
value of the Knudsen number, defined as the ratio of the phonon mean free path
to the characteristic length, $Kn=\frac{\lambda}{L}$.

Figure 1:  Schematic geometry of the studied cases: A thin slab of silicon
device (a) without and (b) with uniform heat generation zone.

##  3 Mathematical Modeling

In 1995, Tzou proposed a non-Fourier approximation for heat conduction, in
which the heat flux vector at a point in the material at time $t+\tau_{\\_}q$
corresponds to the temperature gradient at the same point at time
$t+\tau_{\\_}t$:

$\vec{q}(\vec{r},t+\tau_{\\_}q)=-k\vec{\nabla}T(\vec{r},t+\tau_{\\_}t).$ (9)

Where $\tau_{\\_}q$ and $\tau_{\\_}t$, respectively, represent the heat flux
and the temperature gradient phase lag, both positive and inherent
characteristics of the material. The thermal phase-lag values, like any other
thermal characteristics for engineering materials, can be determined
theoretically or experimentally under various conditions. Tzou [16] has
tabulated some analytical findings for phase lags. Also, Basirat _et al_. [36]
have reported values for the phase lags of some metal films. In the present
research, the DPL model will be called the standard DPL. In 2010, Tzou
generalized the DPL model by introducing the new term measuring the
nonlocality of the heat flux vector. Consequently, the novelly proposed
equation is:

$q(x+\lambda_{\\_}q,t+\tau_{\\_}q)=-K\frac{\partial T}{\partial
x}(x,t+\tau_{\\_}t).$ (10)

The heat flux, $\tau_{\\_}q$, and temperature, ($\tau_{\\_}t$), phase lags,
are considered small in comparison to the time scale, and hence only the first
order Taylor expansion of the DPL equation relative to $\tau_{\\_}q$ and
$\tau_{\\_}t$ are kept. The situation is the same while considering
nonlocality for the heat flux. The correlating length ($\lambda_{\\_}q$) is
also modest compared to the space dimension. Accordingly, the Eq. 10
containing the first-order effects of $\tau_{\\_}q$, $\tau_{\\_}T$, and
$\lambda_{\\_}q$ then yields:

$q(x,t)+\lambda_{\\_}q\frac{\partial q(x,t)}{\partial
x}+\tau_{\\_}q\frac{\partial q(x,t)}{\partial t}=-K\frac{\partial T}{\partial
x}(x,t)-K\tau_{\\_}T\frac{\partial^{2}T}{\partial x\partial t}(x,t).$ (11)

Moreover, the conventional energy equation is,

$-\frac{\partial q}{\partial x}=C\frac{\partial T}{\partial t}.$ (12)

Taking the derivative of the Eq. 11 relative to the position, one can
eliminate the heat flux terms in Eq. 11 using the energy equation,

$(\frac{K}{C\tau_{\\_}q})\frac{\partial^{2}T}{\partial
x^{2}}+(\frac{K\tau_{\\_}T}{C\tau_{\\_}q})\frac{\partial^{3}T}{\partial
x^{2}\partial t}=\frac{1}{\tau_{\\_}q}\frac{\partial T}{\partial
t}+(\frac{\lambda_{\\_}q}{\tau_{\\_}q})\frac{\partial^{2}T}{\partial x\partial
t}+\frac{\partial^{2}T}{\partial t^{2}}.$ (13)

In the next step, to obtain the normalized equations, the following non-
dimensional parameters are defined:

$\theta=\frac{T-T_{\\_}0}{T_{\\_}W-T_{\\_}0},\ t^{*}=\frac{t}{\tau_{\\_}q},\
B=\frac{\tau_{\\_}t}{\tau_{\\_}q},\ \eta=\frac{x}{L},\ Kn=\frac{\lambda}{L},\
\gamma=\frac{\lambda_{\\_}q}{L}.$ (14)

Here, $\lambda$ and $L$, respectively, the phonon mean free path and the
length of the system, determine the Knudsen number, $Kn$. Using the non-
dimensional parameters, the Eq. 13 is rewritten this way:

$\frac{\partial\theta}{\partial
t^{*}}+\gamma\frac{\partial^{2}\theta}{\partial\eta\partial
t^{*}}+\frac{\partial^{2}\theta}{\partial
t^{*^{2}}}=\frac{Kn^{2}}{3}\frac{\partial^{2}\theta}{\partial\eta^{2}}+B\frac{Kn^{2}}{3}\frac{\partial^{3}\theta}{\partial\eta^{2}\partial
t^{*}}.$ (15)

Previously in 2002, Chen [37] simulated the one-dimensional silicon transistor
using the one-dimensional Ballistic-Diffusive Equations and the Boltzmann
equation for the phonons. In the following, the accuracy of our results
obtained utilizing NDPL, will be checked by comparison to Chen’s findings.

###  3.1 Boundary conditions

As the DPL model does not consider the boundary phonon scattering effects,
enforcing no-jump boundary conditions leads to unsatisfactory results,
especially near the boundaries [35]. Hence, here, the temperature jump
boundary condition is used to fix the problem:

$\theta_{\\_}s-\theta_{\\_}w=-\alpha
Kn(\frac{\partial\theta}{\partial\eta})_{\\_}w.$ (16)

Where, $\theta_{\\_}s$ and $\theta_{\\_}w$ are, subsequently, the wall’s jump
temperature and the boundary temperature. Also, $\alpha$, is the coefficient
which will be adjusted to satisfy the boundary condition. Also, two unknown
parameters of $B$ and $\gamma$ alongside $\alpha$, appearing in Eq. 15, are
determined such that the results of the NDPL model becomes compatible with the
solution of the phonon Boltzmann equation (PBE). Interestingly, it is obtained
that the values of these parameters depend on the Knudsen number and time.
Notably, according to the available data, the following values for $\alpha$,
$B$, and $\gamma$ are attained:

$\alpha=\begin{cases}0.35t,&t\leq 0.1\\\ 0.27t,&0.1\leq t\leq 1\\\ 0.55,&t>1\
\&\ Kn>1\\\ 0.65-0.1\times(Kn),&t>1\ \&\ Kn<1.\\\ \end{cases}$ (17)
$B=\begin{cases}0.45t,&t\leq 0.1,\\\ 0.35t,&0.1\leq t\leq 1,\\\
0,&t>1.\end{cases}$ (18) $\gamma=0.35\times(Kn);for\ all\ t.$ (19)

Using the above-obtained constants, results in temperature and heat flux
profiles that are very close to what is calculated from the Boltzmann equation
and also the ballistic-diffusive equations. It should be mentioned that, here,
the thermal properties such as thermal conductivity, sound velocity, and
specific heat are considered constant.

##  4 Numerical solution

The Eq. 13 is solved using a completely implicit first- and second-order
finite difference method. Discretizing of all derivatives is central, where a
stable and convergent three-level finite difference scheme is utilized [38].
Accordingly, the discretizations are written as below:

$\frac{\partial\theta}{\partial t^{*}}=\frac{1}{2\Delta
t}[T^{n+1}_{\\_}i-T^{n-1}_{\\_}i],$ (20)
$\frac{\partial^{2}\theta}{\partial\eta\partial t^{*}}=\frac{1}{4\Delta
x\Delta
t}[(T^{n+1}_{\\_}{i+1}-T^{n+1}_{\\_}{i-1})-(T^{n-1}_{\\_}{i+1}-T^{n-1}_{\\_}{i-1})],$
(21) $\frac{\partial^{2}\theta}{\partial t^{*^{2}}}=\frac{1}{\Delta
t^{2}}[T^{n+1}_{\\_}{i}-2T^{n}_{\\_}{i}+T^{n-1}_{\\_}{i}],$ (22)
$\frac{\partial^{2}\theta}{\partial\eta^{2}}=\frac{1}{\Delta
x^{2}}[T^{n+1}_{\\_}{i+1}-2T^{n+1}_{\\_}i+T^{n+1}_{\\_}{i-1}],$ (23)

and,

$\frac{\partial^{3}\theta}{\partial\eta^{2}\partial t^{*}}=\frac{1}{2\Delta
x^{2}\Delta
t}[(T^{n+1}_{\\_}{i+1}-2T^{n+1}_{\\_}{i}+T^{n+1}_{\\_}{i-1})-(T^{n-1}_{\\_}{i+1}-2T^{n-1}_{\\_}{i}+T^{n-1}_{\\_}{i-1})].$
(24)

Here, a weighted average is implied for stability and convergence. The
convergence of the numerical method is bettered by the Knudsen number
reduction. It is to say that for high Knudsen numbers, the solution explicitly
depends on the marching step size.

##  5 Results and Discussions

The results for modeling the heat transport in 1-D MOSFET using the numerical
simulation of the NDPL model, are manifested in Figs. 2-4. Specifically, the
comparison of the results for the non-dimensional temperature profile and the
heat flux distribution, obtained from the nonlocal DPL, standard DPL, and PBE
are presented. It should be made clear that verification of the output
profiles is performed using the results achieved from PBE. Along the present
NDPL simulation, the profiles for standard DPL modeling are also reproduced by
the authors.

In the work [35], the temperature and heat flux profiles for 1-D MOSFETs with
different Knudsen numbers, at various times, have been presented using the
standard DPL model. While the Knudsen number grows, the deviation of the
results from the precious profiles of the PBE increases. The present study
will try to solve the large Knudsen number issue, using the nonlocal heat
transport concept. More precisely, the nonlocality of heat flux is considered
in the system by introducing the non-dimensional correlating length of
$\gamma$=3.5. This almost large $\gamma$ confirms that taking into
consideration the nonlocality is very crucial for large Knudsen number nano-
structures. Trying to clarify, the better estimate of the results obtained
from the NDPL as compared to that of the standard DPL model is shown in Fig. 2
for the 1-D transistor with the Knudsen number of 10. Specifically, Fig. 2 (a)
presents the temperature profile for the case with Kn=10. As the figure
suggests, the behavior of the temperature distributions obtained from NDPL for
the scaled times of 0.01 and 0.1, are very much like the PBE solutions. It is
obvious that the compatibility of the NDPL-obtained temperature distribution
with PBE results is always better than those attained from DPL. This
consistency is so that the NDPL calculated temperature behavior, completely
fits the PBE result for scaled coordinates below x∗=0.1 and 0.8, respectively,
for t∗=0.01 and 0.1. When t∗=0.01, both the size and time are very low, and
the non-Fourier behavior of the system is very significant. For t∗=0.1, the
system still experiences non-Fourier conditions, but it is a little less in
comparison to that of the t∗=0.01. As evident in Fig. 2 (a), considering
nonlocality in the DPL model, handles non-Fourier behavior in nearly all
positions such that the temperature distribution almost equals the PBE result.
This is true while the result for big times like t∗=10, does not vary
considerably with adding nonlocality characteristics. This department is
justified as although the size of the system is petite, but at significant
times it reaches a steady state, and consequently the non-Fourier behavior is
adjusted.

Figure 2:  The temperature profile calculated from NDPL, DPL, and PBE when
Kn=0.1

The heat flux versus position is demonstrated in Fig. 2 (b). The very better
consistency of the NDPL modeling compared to standard DPL, with the PBE data,
is further seen. When t∗=0.1, for positions below x∗=0.01, similar to the
temperature profile, the heat flux matches well with PBE results. Moreover,
when t∗=0.1 and for x∗¿0.2, the NDPL obtained heat flux profile displays
considerable improvement over what the standard DPL calculates. More
concretely, for large Knudsen numbers, the heat flux distribution obtained
from DPL modeling behaves worse than the computed temperature profile in a way
that gives a much inferior approximation of the exact outcomes of PBE.
Interestingly, adding nonlocality to the DPL model, solves the problem and
leads to the more correct solutions. For larger times, the nonlocality, and
lagging behavior, loses its importance and $\gamma$ and B become zero. This is
expected as the transistors reach the steady state and the role of non-Fourier
heat transfer becomes much less prominent. As Fig. 2 (b) suggests, $\alpha$,
the temperature jump coefficient, is the only responsible parameter for
achieving appropriate temperature and heat flux profiles. In this research,
for large ts, the Knudsen number dependent $\alpha$ is determined such that it
severely rectifies the very different standard DPL obtained heat flux profile.

Figure. 3 shows the state of the 1-D transistor with the Knudsen number of 1.
This Kn is neither high nor low. The trend for temperature and heat flux
behavior is almost similar to that for Kn=10. Here we analyze the results in
more detail. For lower Knudsen numbers, it is anticipated that the non-Fourier
attitude and, consequently the nonlocality effect fades. The non-dimensional
nonlocality presenter, $\gamma$, is 0.35, which is one-tenth the one for
Kn=10. As Fig. 3 (a) confirms, this almost low $\gamma$, is essential in
getting results that are more congruent to the PBE solutions. Also, as non-
Fourier behavior has become a little less important in comparison to the case
with Kn=10, considering $\gamma$ makes the NDPL calculated profiles to be very
accurately fitted to the available precise data [37]. In other words, almost
the complete deviation of the standard DPL results from PBE outputs’, is
filled via the nonlocality contribution. For instance, when t∗=0.1, for a
broader range of position, say for x∗¡1.5 and x∗¿2, the NDPL temperature
profile finding, provides the results found out solving the PBE. The
consistency is more noticeable when we are dealing with t∗=1. For this time,
compatibility and fitness exist for almost all positions. In consequence,
taking into account the temperature profile obtained from standard DPL, it
becomes pretty clear how our new model is doing much better. Although, for
Kn=10, the NDPL distributions were close to the PBE solutions, the better
agreeableness at Kn=1.0, is attributed to the slightly weaker non-Fourier
behavior, which is practically covered by the simultaneous consideration of
nonlocality and phase lagging. When a long time passes, the situation becomes
like in nano-structures with larger Kn numbers. In similarity, the parameters
$\gamma$ and $B$ reach zero, and the effect of nonlocality and phase lagging
disappears.

On the other side, in Fig. 3 (b), the heat flux demeanor versus position is
illustrated. Same as the case with Kn=10, here also, the NDPL modeling assures
the results very close to the available atomistic data. As it is viewed, when
t∗=0.1, like what happens to temperature profile, the NDPL calculated heat
flux, is consistent with available data for x∗¡0.1. Only as it is the case for
Kn=10, the heat flux at x∗=0, is a bit larger than that of PBE. Moreover, heat
flux fits preciously for x∗¿0.2, and well for x∗¡0.2, with the existing data.
One can conclude that although NDPL predicts the heat flux correctly, a slight
deviation exists for meager non-dimensional positions. It is important to note
that the accuracy of the heat flux profile obtained from standard DPL is
reported to be always less than the DPL calculated temperature distribution,
while the accuracy is defined in closeness to the PBE results. Intriguingly,
the NDPL resolves the heat flux inaccuracy and gives a highly meticulous
distribution. Also, for a larger scaled time of t∗=10, when $\gamma$ and B are
zero, the newly defined Kn dependent $\alpha$, reproduces the PBE solution.

Figure 3:  The same as Fig. 2 but for Kn=1.

At last, the temperature and heat flux distributions are studied for low
Knudsen numbers. The results are shown in Fig. 4. As low Knudsen number
systems have larger sizes, one anticipates less anomalous heat transfer
inside. This does not mean that the effect of nonlocality on heat flux
disappears. Even though the non-dimensional nonlocality of the heat flux for
Kn=0.1, is $\gamma$=0.035, which is two orders of magnitude less than
$\gamma$=3.5 for Kn=10, it is crucial to consider such nonlocality to get more
precise plots. As Fig. 4(a) suggests for the lower scaled time of t∗=1, where
the Non-Fourier attitude is bold compared with larger times of t$*$=10 and
t∗=100, the NDPL presents a temperature profile with signifiently much better
consistency with the PBE relative to the standard DPL. On the other hand, the
small Knudsen numbers and large scaled times, say t∗¿1, change the conditions
in favor of removing the nonlocality and phase lag effects. So, while $\gamma$
and B become zero, the MOSFET reaches the steady state, and consequently, both
the DPL and NDPL calculate distributions nearly the same as the PBE findings.

The imprecision in calculating the heat flux for low Kn number and smaller
times using the DPL model is also worked out, implying the newly proposed
NDPL. Just as it is evident in Fig. 4, for kn=0.1 and t∗=1, the NDPL obtained
heat flux presents a remarkably improved condition such that at x∗=0 is closer
to the PBE result and also fits it for a much wider range of positions. For
larger times, as previously mentioned for temperature distribution, both DPL
and NDPL work well while the non-Fourier effect is adjusted.

Figure 4:  Comparison of the temperature distribution obtained from the
Boltzmann Equation, right triangle, the nonlocal behavior DPL model, dashed
line, and the standard DPL model for different instantaneous times and Kn=0.1.

##  6 Conclusion

The current research tries to establish a framework for micro/nanoscale heat
transport investigation that produces precious results while concurrently
having low computational costs. To achieve this aim, the 1-D nonlocal DPL
model with the possibility of generalization to the higher dimensions is
developed. This method considers the nonlocality of the heat flux along the
phase lagging, while the temperature jump is also enforced. In particular, a
non-dimensional parameter $\gamma$, in addition to the $\alpha$, the
temperature jump coefficient, and $B$, the phase lagging ratio, both defined
in DPL, is introduced to determine the strength of the nonlocality. Then, the
formulated computational procedure is applied to a 1-D silicon MOSFET. The
parameters $\gamma$, $\alpha$, and $B$ are newly evaluated by verifying of the
obtained results with PBE-calculated available data. The presence of the
Knudsen-dependent nonlocality coefficient, $\gamma$, is found to be crucial
for obtaining accurate temperature and heat flux distribution. In further
detail, as expected, for the large Knudsen number, the nonlocality at heat
flux is very distinguished. As an instance, when Kn is 10, the $\gamma$
parameter becomes 3.5, which is one hundred times the $\gamma$ value for
Kn=0.1. This does not mean that nonlocality is not essential for low Knudsen
numbers. Taking into account the small $\gamma$ also makes the results very
close to the PBE available data, especially for earlier times. In addition,
introducing nonlocality to the DPL model, resolves the problem of inaccurate
heat flux obtained from DPL and corrects it to remarkably acceptable value.
Finally, when the heat transport inside a MOSFET is predicted well, it is more
feasible to propose new transistors supporting the thermal management
solutions.

## References

  * [1] R. Mahajan, R. Nair, V. Wakharkar, J. Swan, J. Tang, G. Vandentop, Emerging directions for packaging technologies, Intel Technol. J. 6 (2002) 2.
  * [2] M. Moghaddam, J. Ghazanfarian, A. Abbassi, Implementation of DPL-DD model for the simulation of nanoscale MOS devices, IEEE Transactions on Electron Devices, 61(9), 3131-3138, 2014.
  * [3] A.J. Minnich, Determining phonon mean free paths from observations of quasiballistic thermal transport, Physical review letters, 109(20), 205901, 2012.
  * [4] J. Ghazanfarian, Z. Shomali, S. Xiong, 21st Century Nanoscience–A Handbook: Nanophysics Sourcebook (Volume One), Sattler, K. D. (Ed.), Chapter 4, CRC Press, 2019.
  * [5] Y. H. Chiu, V. V. Deshpande, H. C. Postma, C. N. Lau, C. Miko, L. Forro, and M. Bockrath, Ballistic phonon thermal transport in multiwalled carbon nanotubes, Physical review letters, 95(22), 226101, 2005.
  * [6] F. X. Alvarez, and D. Jou, Memory and nonlocal effects in heat transport: From diffusive to ballistic regimes, Applied physics letters, 90(8), 083109, 2007.
  * [7] Z. Shomali, R. Asgari, Effects of low-dimensional material channels on energy consumption of nano-devices, International Communications in Heat and Mass Transfer, 94, 77-84, 2018.
  * [8] T. Q. Qiu and C. L. Tien,Short-pulse laser heating on metals, International Journal of Heat and Mass Transfer, 35, 719, 1992.
  * [9] R. A. Guyer and J. A. Krumhansl, Solution of the linearized Boltzmann equation, Physical Review, 148, 766, 1966.
  * [10] M. E. Gurtin and A. G. Pipkin, A general theory of heat conduction with finite wave speed,  Archive for Rational Mechanics and Analysis, 31,113, 1968.
  * [11] W. Roetzel, N. Putra, S. K. Das,  International Journal of Thermal Science 42 (2003) 541.
  * [12] H. Basirat Tabrizi, S. Andarwa,  International Communication in Heat and Mass Transfer 36 (2009) 186.
  * [13] G. D. Mahan, and F. Claro, Nonlocal theory of thermal conductivity. Physical Review B 38(3) (1988) 1963.
  * [14] B. Vermeersch, and A. Shakouri, Nonlocality in microscale heat conduction, arXiv preprint arXiv:1412.6555 (2014).
  * [15] Da Yu Tzou,  The generalized lagging response in small-scale and high-rate heating,  International Journal of Heat and Mass Transfer 38(17) (1995) 3231.
  * [16] D.Y. Tzou, Macro- to Microscale Heat Transfer: The Lagging Behavior. Taylor & Francis, Washington, D.C., USA, 1997.
  * [17] D. Y. Tzou, and Z. Y. Guo, Nonlocal behavior in thermal lagging,  International Journal of Thermal Sciences 49(7) (2010) 1133.
  * [18] B. Y. Cao, Z. Y. Guo, Equation of motion of a phonon gas and non-Fourier heat conduction,  J. Appl. Phys. 102 (2007) 053503.
  * [19] J. Ghazanfarian, Z. Shomali, A. Abbassi, Macro-to nanoscale heat and mass transfer: the lagging behavior, Int. J. Thermophys 36 (2015) 1416.
  * [20] Z. Shomali, R. Kovács, P. Ván, I.V. Kudinov, and J. Ghazanfarian,  Recent Progresses and Future Directions of Lagging Heat Models in Thermodynamics and Bioheat Transfer,  Continuum Mechanics and Thermodynamics, 34:637–679, 2022.
  * [21] R.S. Samian, A. Abbassi, J. Ghazanfarian, Thermal investigation of common 2d FETs and new generation of 3-d FETs using Boltzmann transport equation in nanoscale, Int. J. Mod. Phys. C 24 (2013) 1350064.
  * [22] R.S. Samian, A. Abbassi, J. Ghazanfarian, Transient conduction simulation of a nanoscale hotspot using finite volume lattice Boltzmann method, Int. J. Mod. Phys. C 25 (04) (2014) 1350103.
  * [23] A.L. Moore, L. Shi, Emerging challenges and materials for thermal management of electronics, Mater. Today 17 (4) (2014) 163.
  * [24] J. Ghazanfarian and Z. Shomali,  Investigation of dual-phase-lag heat conduction model in a nanoscale metal-oxide-semiconductor field-effect transistor,  International Journal of Heat and Mass Transfer, 55(21-22):6231-7, 2012.
  * [25] Z. Shomali, J. Ghazanfarian, A. Abbassi,  Investigation of bulk/film temperature-dependent properties for highly non-linear DPL model in a nanoscale device: the case with high-k metal gate MOSFET,  Superlattices and Microstructures, 83:699, 2015.
  * [26] Z. Shomali, A. Abbassi, J. Ghazanfarian, Development of non-Fourier thermal attitude for three-dimensional and graphene-based MOS devices, Appl. Therm. Eng. 104 (2016) 616.
  * [27] Z. Shomali, B. Pedar, J. Ghazanfarian, A. Abbassi,  Monte-Carlo Parallel Simulation of Phonon Transport for 3D Nano-Devices,  International Journal of Thermal Sciences, 114:139-154, 2017.
  * [28] Z. Shomali, J. Ghazanfarian, A. Abbassi, 3-D Atomistic Investigation of Silicon MOSFETs, In Proceedings of CHT-17 ICHMT International Symposium on Advances in Computational Heat Transfer, ICHMT Digital Library Online, Begel House Inc., 2017.
  * [29] M. H. Fotovvat and Z. Shomali,  A time-fractional dual-phase-lag framework to investigate transistors with TMTC channels (TiS${}_{\\_}3$, In${}_{\\_}4$Se${}_{\\_}3$) and size-dependent properties, 168, 207304, 2022.
  * [30] Z. Shomali, An investigation into the reliability of newly proposed MoSi${}_{\\_}2$N${}_{\\_}4$/WSi${}_{\\_}2$N${}_{\\_}4$ field effect transistor: A monte carlo study, arXiv preprint arXiv:2305.04327 (2023).
  * [31] E. Pop, R.W. Dutton, K.E. Goodson, Monte Carlo simulation of joule heating in bulk and strained silicon, Appl. Phys. Lett. 86 (2005) 082101.
  * [32] E. Pop, S. Sinha, K.E. Goodson, Heat generation and transport in nanometer-scale transistors, Proc. IEEE 94 (8) (2006) 1587.
  * [33] E. Pop, Energy dissipation and transport in nanoscale devices, Nano Res. 3 (3), (2010) 147.
  * [34] S. Gong, L. Chen, H. Feng, Z. Xie, F. Sun, Constructal optimization of cylindrical heat sources surrounded with a fin based on minimization of hot spot temperature, Int. Commun. Heat Mass 68 (2015) 1.
  * [35] J. Ghazanfarian, A. Abbassi,  Effect of boundary phonon scattering on Dual-Phase-Lag model to simulate micro-and nanoscale heat conduction,  International Journal of Heat and Mass Transfer, 52(15-16) (2009) 3706.
  * [36] H. Basirat, J. Ghazanfarian, P. Forooghi, Implementation of dual-phase- lag model at different Knudsen numbers within slab heat transfer, in: Proceedings of International Conference on Modeling and Simulation (MS06), August, 2006, Konia, Turkey, pp. 895899.
  * [37] G. Chen,  Ballistic-diffusive equations for transient heat conduction from nano to macroscale,  ASME J. Heat Transfer 124 (2001) 320–328.
  * [38] W. Dai, L. Shen, R. Nassar, and T. Zhu, A stable and convergent three-level finite difference scheme for solving a dual-phase-lagging heat transport equation in spherical coordinates,  Int. J. Heat Mass Transfer 47 (2004) 1817.


11institutetext: Laboratoire de Physique de l’Ecole Normale Supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris,
France
11email<EMAIL_ADDRESS>22institutetext: INAF–Osservatorio
Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
33institutetext: Istituto di Astrofisica e Planetologia Spaziali (IAPS), INAF,
Via Fosso del Cavaliere 100, 00133 Roma, Italy

# The Orion-Taurus ridge: a synchrotron radio loop at the edge of the Orion-
Eridanus superbubble

Andrea Bracco 11 Marco Padovani 22 and Juan D. Soler 33

(Received June 26, 2023; accepted August 31, 2023)

Large-scale synchrotron loops are recognized as the main source of diffuse
radio-continuum emission in the Galaxy at intermediate and high Galactic
latitudes. Their origin, however, remains rather unexplained.

Using a combination of multi-frequency data in the radio band of total and
polarized intensities, for the first time in this letter, we associate one arc
– hereafter, the Orion-Taurus ridge – with the wall of the most prominent
stellar-feedback blown shell in the Solar neighborhood, namely the Orion-
Eridanus superbubble.

We traced the Orion-Taurus ridge using 3D maps of interstellar dust extinction
and column-density maps of molecular gas, $N_{\rm H_{2}}$. We found the Orion-
Taurus ridge at a distance of 400 pc, with a plane-of-the-sky extent of $180$
pc. Its median $N_{\rm H_{2}}$ value is $(1.4^{+2.6}_{-0.6})\times 10^{21}$
cm-2.

Thanks to the broadband observations below 100 MHz of the Long Wavelength
Array, we also computed the low-frequency spectral-index map of synchrotron
emissivity, $\beta$, in the Orion-Taurus ridge. We found a flat distribution
of $\beta$ with a median value of $-2.24^{+0.03}_{-0.02}$ that we interpreted
in terms of depletion of low-energy ($<$ GeV) cosmic-ray electrons in recent
supernova remnants ($10^{5}$ \- $10^{6}$ yrs). Our results are consistent with
plane-of-the-sky magnetic-field strengths in the Orion-Taurus ridge larger
than a few tens of $\mu$G ($>30-40\,\mu$G).

We report the first detection of diffuse synchrotron emission from cold-
neutral, partly molecular, gas in the surroundings of the Orion-Eridanus
superbubble. This observation opens a new perspective to study the multiphase
and magnetized interstellar medium with the advent of future high-sensitivity
radio facilities, such as the C-Band All-Sky Survey and the Square Kilometre
Array.

###### Key Words.:

magnetic fields – radiation mechanisms: non-thermal – ISM: cosmic rays – ISM:
clouds – ISM: dust, extinction – radio continuum: ISM

##  1 Introduction

The intermediate-to-high Galactic latitude sky brightness in radio continuum
is dominated by numerous large-scale loops that have long been observed (e.g.,
Berkhuijsen, 1971; Haslam et al., 1971) and interpreted as synchrotron
emission structures (Vidal et al., 2015; Planck 2015 results. XXV., 2016;
Thomson et al., 2021; Martire et al., 2023), but whose origin remains
uncertain. Most of them have been associated with remnants of supernova (SN)
explosions (Berkhuijsen et al., 1971; Salter, 1983), where cosmic-ray
electrons (CRe) are accelerated in shock-compressed magnetic fields in the
interstellar medium (ISM), giving rise to synchrotron emission (Ferriere &
Zweibel, 1991). The superposition of such loops from SN remnants may also be
responsible for most of the total radio emission along the Galactic plane
itself, as postulated by Mertsch & Sarkar (2013). However, to date, a
reference model on the physical origin of synchrotron radio loops is still
lacking (e.g., Sofue, 2015; Kataoka et al., 2021; Panopoulou et al., 2021).

In this letter, we present the emission and the physical properties of the
Orion-Taurus ridge, which is a 25-degree-wide radio loop located in the
constellation of Orion and centered around Galactic coordinates
$(l,b)=(190^{\circ},-20^{\circ})$. Despite being already detected (e.g.,
Milogradov-Turin & Urošević, 1997; Borka et al., 2008; Vidal et al., 2015;
Planck 2015 results. XXV., 2016), the Orion-Taurus ridge represents a rather
unreported radio structure compared to more renown loops, such as Loop I, II,
III, and IV (Berkhuijsen, 1971; Vidal et al., 2015).

For the first time, we show that the Orion-Taurus ridge is a synchrotron loop
associated with cold neutral ISM at the northern edge of the Orion-Eridanus
superbubble (OES, Heiles, 1976; Reynolds & Ogden, 1979; Ochsendorf et al.,
2015). The OES is one of the closest-by (distance, $D<500$ pc) and largest
circular structures in the sky (see the H$\alpha$ map from Finkbeiner (2003)
in the right panel of Fig. 1). It includes the Orion molecular cloud (OMC) and
is considered the joint result of strong ionizing UV radiation, stellar winds,
and SN explosions from the Orion OB1 association (e.g., Heiles, 1997; Bally,
2008).

Figure 1:  Stereographic projections showing the Orion-Taurus ridge. Left:
Brightness temperature ($T_{\rm b}$) at 74 MHz observed by the LWA telescope
at 2∘ FWHM resolution. The contours are at 2100, 2200, 2400, 2600, and 2800 K,
respectively. Center: Polarized intensity ($PI$) at 30 GHz observed by the
Planck satellite at the same FWHM resolution of $T_{\rm b}$(74 MHz) with the
plane-of-the-sky magnetic field orientation traced by line integral
convolution (Cabral & Leedom, 1993). Black contours are the same as in the
left panel. Right: H$\alpha$ emission map from Finkbeiner (2003) showing the
Orion-Taurus ridge in contours and the extent of the Orion-Eridanus
superbubble (OES) traced in blue as in Soler et al. (2018). The same Galactic
coordinate grid is shown in all panels with steps in $l$ and $b$ of
$10^{\circ}$ and centered in ($l,b$) = ($200^{\circ},-30^{\circ}$). The
stereographic projection was produced with the Python routine
stereo$\\_$proj$\\_$hpx that can be found at
http://github.com/abracco/cosmicodes/blob/master/4GMIMS/Planck_routines.py.

In the following, we introduce the multi-frequency data used in the analysis
(Sect. 2). We detail the morphological correlation between radio data of total
and polarized intensities tracing synchrotron in the Orion-Taurus ridge and
the dusty, partly molecular, ISM at the edge of the OES (Sect. 3). We compute
the spectral index of synchrotron emissivity, $\beta$, in the Orion-Taurus
ridge using radio continuum data below 100 MHz (Sect. 4). Based on the
parametric model of Padovani et al. (2018) for the energy spectrum of CRe
responsible for the low-frequency radio emission (Padovani et al., 2021), we
propose an original mechanism to explain the values of $\beta$ in the Orion-
Taurus ridge.

##  2 Data

In this section, we describe the multi-wavelength data set used to probe the
Orion-Taurus ridge.

###  2.1 Synchrotron total emission: Long Wavelength Array data

We used data from the northern-sky survey of the Long Wavelength Array
(LWA)111http://lda10g.alliance.unm.edu/LWA1LowFrequencySkySurvey/ as described
in Dowell et al. (2017, hereafter D17). The original brightness
temperature,$T_{\rm b}(\nu)$, maps and their uncertainties correspond to nine
frequency bands between 35 and 80 MHz with angular resolution between FWHM
$\approx 5^{\circ}$ and $\approx 2^{\circ}$, respectively (see Table 1 in
D17). These maps trace diffuse synchrotron emission in the Galaxy.

In this work, we only focused on the five highest frequency bands, namely at
50, 60, 70, 74, and 80 MHz, in order to reach FWHM of $w_{0}=3.5^{\circ}$, at
worst. All LWA maps, in units of K, are downloadable in
HEALPix222http://healpix.sf.net format (Górski et al., 2005) with
$N_{\textrm{s}ide}$ $=$ 256 (13.7$\arcmin$ wide pixels). In Sect. 4, where we
study the synchrotron spectral index of the loop, we smoothed the maps to the
common resolution of $w_{0}$, assuming Gaussian beams333The effective FWHM is
$\sqrt{w_{0}^{2}-w_{i}^{2}}$, where $w_{i}$ represents the Gaussian FWHM of
each original map. and using the healpy.sphtfunc.smoothing Python package
(Zonca et al., 2019).

The Orion-Taurus ridge can be seen in the left panel of Fig. 1, where, as an
example, we show it with white contours on top of the $T_{\rm
b}(74\leavevmode\nobreak\ {\rm MHz})$ map. The Orion-Taurus ridge was isolated
by unsharp masking the $T_{\rm b}$(74 MHz) map using a Gaussian window of FWHM
$=$ $10^{\circ}$ (see also Vidal et al., 2015), and masking out all filtered
values below 40 K. The choice of this threshold was only used to optimize the
sky selection of the Orion-Taurus ridge that was finalized with the function
label of the scipy Python package. The Orion-Taurus ridge extends over several
tens of degrees both in Galactic longitude and latitude.

###  2.2 Synchrotron polarization: Planck data

Because of the numerous instrumental and physical challenges that affect low-
frequency radio-polarimetric data, such as those of LWA (Taylor et al., 2012)
– e.g., Stokes-parameter leakage, ionospheric and interstellar Faraday
rotation – it was easier to spot the Orion-Taurus ridge in synchrotron
polarization at high frequency from satellite observations. We inspected the
30 GHz Planck maps of Stokes $Q$ and $U$ (Planck 2018 results. I,
2020)444http://pla.esac.esa.int/ for the Orion-Taurus ridge, after smoothing
them at FWHM $=$ $2^{\circ}$. In the right panel of Fig. 1, we show the Orion-
Taurus ridge with black contours on the corresponding polarized intensity
($PI$ $=$ $\sqrt{Q^{2}+U^{2}}$). It was already discovered by Planck 2015
results. XXV. (2016) but never discussed in depth. Although $PI$ extends more
toward the east, the visual correlation with most of the loop seen in $T_{\rm
b}$(74 MHz) supports its origin as a synchrotron structure. Like most of the
radio loops (Vidal et al., 2015), the Orion-Taurus ridge also appears
elongated along the plane-of-the-sky (POS) magnetic field, $\@vec{B}_{\perp}$,
which we inferred from synchrotron polarization rotating by 90∘ the 30-GHz
polarization angle, $\psi$ $=$ $0.5\arctan(-U,Q)$ 555The minus sign converts
the Planck Stokes U to the IAU convention.. The $\@vec{B}_{\perp}$ orientation
is shown using the line integral convolution (LIC, Cabral & Leedom, 1993) from
healpy.

In this work, we considered Planck data to trace the total gas column density
($N_{\rm H}$) along the line of sight (LOS). In particular, we used the
$N_{\rm H}$ map derived from the dust optical depth at 353 GHz ($\tau_{353}$),
which was obtained by fitting a modified blackbody spectrum to the Planck
intensity data below 857 GHz and IRAS at 100 $\mu$m (Planck 2013 results. XI.,
2014). The $\tau_{353}$ map (FWHM = 5$\arcmin$) was converted into $N_{\rm H}$
using a conversion factor of $8.3\times 10^{25}$ cm-2 (Planck Collaboration
Int. XXXV, 2016). The $N_{\rm H}$ map is shown in a composite image in Fig. 2.

Finally, as discussed in Sect. 4, we used the Planck emission measure ($EM$)
map to estimate the contribution of free-free absorption to synchrotron
emission. The maximum-likelihood (ML) value of the $EM$
map${}^{\rm\ref{planckdata}}$ was the result of a Bayesian analysis for
component separation (Planck 2015 results. X., 2016). We smoothed the original
map to the resolution of $w_{0}$.

###  2.3 Orion-Taurus ridge distance: 3D dust map

We made use of the publicly available 3D dust map from Lallement et al.
(2019)666http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/625/A135 to constrain
the distance to the Orion-Taurus ridge. This data cube was obtained by
combining _Gaia_ DR2 and 2MASS photometric data with _Gaia_ DR2 parallaxes.
The LOS dust distribution was inferred with a regularized Bayesian tomographic
inversion of starlight extinction measurements. This technique has the
advantage of accounting for spatial correlation of dust structures with a
limited physical resolution of 25 pc within 1 kpc from the Sun. For our
purpose, however, we did not require high resolution in the 3D dust
reconstruction (e.g., Green et al., 2019; Leike et al., 2020). The 25 pc
spacing corresponds to the value of $w_{0}$ at a distance of 400 pc, which is
approximately the distance of the OMC (e.g., Großschedl et al., 2021). The 3D
dust map contains values of differential extinction (${\rm d}A_{\rm V}/{\rm
d}D$) at 5500 Å in units of ${\rm mag\leavevmode\nobreak\ pc^{-1}}$. Examples
of the integrated 3D dust maps are shown in Fig. 2 and Fig. 6.

##  3 Morphological correlation analysis

Figure 2:  Stereographic view of the 3D dust map integrated between 300 and
500 pc (in blue) and the total gas column density ($N_{\rm H}$, in green)
toward the Orion-Taurus ridge, shown with white contours. Light-blue contours
trace the 3D dust map at $A_{\rm V}=2,3,4,5$ mag. The $N_{\rm H}$ map is shown
above $10^{21}$ cm-2. The coordinate grid and the white contours are the same
as in Fig. 1. The long-dashed yellow line indicates the edge of the OES. The
short-dashed-yellow line indicates the edge of the Orion shell reported by
Foley et al. (2023). A white circle highlights the center of the OES.

The Orion-Taurus ridge is seen both in $T_{\rm b}$(74 MHz) and in $PI$(30 GHz)
(see Fig. 1). It is a synchrotron structure with curvature directed toward the
center of the OES, or possibly of a subset of it at a comparable distance
named the Orion shell (see Fig. 2 and Foley et al., 2023).

In Fig. 2 we show a composite image, where the tentative edge of the
superbubble (e.g., Soler et al., 2018) is sketched in yellow on top of both
the $N_{\rm H}$ map (in green) and the 3D dust map integrated between 300 and
500 pc (${A}_{\rm V}(300,500{\rm pc})$, in blue). The Orion-Taurus ridge,
delimited by white contours as in Fig. 1, closely follows the northern dusty
edge of the OES. It reveals a coherent morphological correlation with the 3D
dust map extending from the OMC (A and B) to the more diffuse western rim with
${A}_{\rm V}(300,500{\rm pc})\gtrsim 5$ mag and $N_{\rm H}\gtrsim 10^{21}$
cm-2. The mean value of ${A}_{\rm V}(300,500{\rm pc})$ within the Orion-Taurus
ridge is 8 mag, which integrated over 200 pc corresponds to a lower
limit777This is only a lower limit because most likely the gas distribution is
not uniform within 200 pc, but condensed on a smaller length scale. to the
mean number density ($n_{\rm H}$) of $n_{\rm H}\approx 50$ cm-3, typical of
the cold neutral medium (Ferrière, 2020). As shown in Appendix A, this visual
correlation would not be attained if the 3D dust map was integrated in
different distance ranges, between 100 and 300 pc or between 500 and 700 pc
(see Fig. 6).

We quantified the correlation by comparing the radial profiles of the 3D dust
map with that of $T_{\rm b}$(74 MHz) from the center of the OES (see the white
circle in Fig. 2). We defined the plane-of-the-sky distance ($r^{\rm MAX}$, in
units of degrees) from the center of the superbubble to the first peak of a
given radial profile. The peak was found using the find-peaks function of
scipy. Note that $r^{\rm MAX}$ was only considered if the peak of the radial
profiles was at least 20% of the maximum value of the corresponding map. In
Fig. 3, we show a polar plot encoding the distributions of $r^{\rm MAX}$ as a
function of the azimuthal angle increasing from $0^{\circ}$ to $360^{\circ}$
following the E-W direction. We binned the azimuthal angle in steps of
$5^{\circ}$. The plot compares $r^{\rm MAX}$ defined at 74 MHz (orange) with
that of the 3D dust map (blue). The Orion-Taurus ridge only correlates, over
$135^{\circ}$ between $45^{\circ}$ and $180^{\circ}$, with the structure of
the 3D dust map toward the northern rim of the OES at about 400 pc (see
Appendix A for comparison with other distance ranges). For this particular
case, in the right panel of Fig. 3, we also show the difference between the
two $r^{\rm MAX}$ distributions, which peaks at $0^{\circ}$ with a dispersion
smaller than $5^{\circ}$. The median value, the 16th, and 84th percentile
errors are $\Delta r^{\rm MAX}_{\rm m}=0.0^{+1}_{-4}$ ∘.

Our morphological analysis established the correlation between the Orion-
Taurus ridge, seen in total and polarized synchrotron intensities, with cold
neutral medium traced by dust extinction. We also checked the molecular gas
column density map, $N_{{\rm H}_{2}}$, inferred from observations of atomic
hydrogen, HI, emission spectra presented in Kalberla et al. (2020). The Orion-
Taurus ridge shows up at $N_{{\rm H}_{2}}\gtrsim 10^{21}$ cm-2 with a median
value of $(1.4^{+2.6}_{-0.6})$ $\times$ $10^{21}$ cm-2. To our knowledge, this
is the first detection of diffuse synchrotron emission from molecular gas
along the northern edge of the OES.

Figure 3: Morphological correlation between the Orion-Taurus ridge and the 3D
dust map. Left panel: Polar plot of $r^{\rm MAX}$ as a function of the
azimuthal angle in bins of 5∘ (increasing in the E-W direction). $r^{\rm MAX}$
is computed for $T_{\rm b}$(74 MHz) (orange) and the 3D dust map integrated
between 300 and 500 pc (blue). Right panel: Histogram of the difference
between the $r^{\rm MAX}$ distributions shown on the left. The reported values
correspond to the median and the corresponding errors (16th and 84th
percentiles).

##  4 Synchrotron spectral index map

Thanks to the broadband frequency coverage of the LWA, we could study in
detail the synchrotron properties of the Orion-Taurus ridge. We used low-
frequency data as synchrotron represents the only relevant mechanism in this
radio range, without the free-free contamination in emission typical of the
high-frequency end ($>$ 1 GHz). Given the five LWA bands introduced in Sect.
2.1, we computed the spectral index of synchrotron emissivity, $\beta$ $=$
${\rm d}\ln(T_{\rm b}(\nu))/{\rm d}\ln(\nu)$. The map is shown in Fig. 4,
where all the displayed values have a signal-to-noise ratio larger than 5. Our
$\beta$ map is consistent with that reported in D17 at a lower angular
resolution, which was obtained based on all nine LWA bands (see figure 11).
These values of $\beta$ are slightly flatter than those found at 5∘ FWHM by
Guzmán et al. (2011) – ranging between $-2.35$ and $-2.45$, see their Fig. 6.
Despite this difference, as in Guzmán et al. (2011), we also found a similar
gradient between the main body of the arc and the location of the OMC. As
shown by the blue histogram on the top of Fig. 4, the OMC peaks at
$\beta=-2.15$ with generally flatter values than the rest of the arc, which
has a median value of $\beta_{\rm m}=-2.21^{+0.06}_{-0.03}$.

Figure 4: Synchrotron spectral-index map ($\beta$, FWHM $\approx 3.5^{\circ}$,
colors) of the Orion-Taurus ridge computed using LWA data at 50, 60, 70, 74,
and 80 MHz. The $\beta$ map is overlaid on the $N_{\rm H}$ map from Planck.
The inset above shows the histograms of $\beta$ before (blue) and after
(orange) free-free absorption correction. The OMC corresponds to the blue peak
at $\beta=-2.15$.

###  4.1 $\beta$-flattening from free-free

One possible source of flattening of the synchrotron spectral index below 100
MHz comes from contamination of bremsstrahlung, free-free, absorption along
the LOS (Dowell et al., 2017; Stanislavsky et al., 2023). This effect was
first hypothesized in studies of radio spectra of SN remnants (Kassim, 1989)
and proposed to infer the properties of the ionized gas in the Galaxy from
extended envelopes of H ii regions to the warm ionized medium (WIM, Reynolds
et al., 1998). Because of the scattering of radio photons with ionized gas
along the LOS the observed synchrotron emission may be reduced by a factor
$e^{-\tau_{\rm ff}(\nu)}$, where $\tau_{\rm ff}(\nu)$ represents the free-free
optical depth at frequency $\nu$. Following Stanislavsky et al. (2023), this
can be expressed as

$\tau_{\rm ff}(\nu)=3.014\times 10^{4}Z\left(\frac{T_{e}}{[{\rm
K}]}\right)^{-1.5}\left(\frac{\nu}{[{\rm
MHz}]}\right)^{-2}\left(\frac{EM}{[{\rm cm^{-6}\,pc}]}\right)g_{\rm ff},$ (1)

where $Z$ is the average number of ion charges, $T_{e}$ is the electron
temperature, and $EM$ is the emission measure, that is, the integral along the
LOS of the electron density, $n_{e}$, squared (i.e., $EM=\int_{\rm
LOS}n^{2}_{e}{\rm d}l$). The $g_{\rm ff}$ value is the Gaunt factor that we
assume equal to $\ln{(49.55/Z/\nu)}+1.5\ln{T_{e}}$ (Stanislavsky et al.,
2023).

Considering the vicinity of the Orion-Taurus ridge with the Orion OB1
association, where most of the hydrogen ($Z=1$) must be photodissociated by
strong UV radiation, we applied Eq. 1 with an average $T_{e}=10^{4}$ K
(Ochsendorf et al., 2015) and the $EM$ map from $\it Planck$ (see Sect. 2).
After correction of the LWA $T_{\rm b}$ maps, we derived the spectral index
distribution shown in orange in Fig. 4, with a median value of $\beta^{\rm
c}_{\rm m}=-2.24^{+0.03}_{-0.02}$, where ”c” stands for “corrected”. If the
OMC peak is strongly suppressed after the free-free correction, the average
$\beta$ distribution remains the same. We propose another mechanism
contributing to lowering the observed $\beta$ in LWA data, which relates to
variations of the CRe energy spectrum possibly caused by interacting shocks in
the Orion-Taurus ridge and variability of emission sources.

###  4.2 Proposed mechanism for explaining $\beta$

At a given frequency, the synchrotron spectral index depends on both the CRe
energy spectrum and the strength of $\@vec{B}_{\perp}$. Synchrotron emission
below a few hundred MHz is mostly produced by CRe with energies ¡ GeV for
$B_{\perp}$ of a few $\mu$G (Padovani et al., 2021). Their exact energy
spectrum, however, remains uncertain. Considering local measurements of the
CRe energy flux (e.g., Orlando, 2018), and typical $B_{\perp}$ in the ISM, it
is difficult to account for the flat $\beta$ found in this work (see Fig. 9 in
Padovani et al., 2021, where $\beta\leq-2.48$ at 60 MHz for $B_{\perp}$
$\approx$ 30 $\mu$G).

We modelled the low-energy part of the CRe spectrum using the parametric
approach presented in Padovani et al. (2018) and fixed the high-energy part (¿
10 GeV) to the interstellar value. This choice is supported by high-energy
measurements of the CR-proton spectrum in the OES with the Fermi-LAT satellite
that found no significant difference with respect to the diffuse ISM (Joubaud
et al., 2020).

At low energy, as detailed in Appendix B, using Eq. 2 we could depart from
standard interstellar conditions (Orlando, 2018) and control the population of
CRe. We were only able to reproduce the range of $\beta$ observed with LWA by
severely damping the population of low-energy CRe ($\sim$100 MeV to GeV) and
considering strong $B_{\perp}$ larger than a few tens of $\mu$G. These values
of $B_{\perp}$ are almost one order of magnitude larger than the average LOS
magnetic-field strength of $\approx$ $-3\,\mu$G measured with Zeeman splitting
of HI lines toward the OMC (Heiles, 1997).

Although these results highly depend on our limited knowledge of the CRe
energy spectrum in the Galaxy, the depletion of low-energy CRe can be related
to the specific conditions of the Orion-Taurus ridge as part of the OES (see
Sect. 3). In the context of CRe properties in SN remnants, recent work by
Morlino & Celli (2021) showed that a spectral break in the CRe energy spectrum
below a few GeV should be expected as a consequence of the drop in
acceleration efficiency of the CRe during the adiabatic stage of the remnant
expansion. This is a plausible explanation as the OES is considered the result
of highly episodic star formation in the last 10-15 Myrs. Several SN
explosions contributed to shaping the cavity and driving the global dynamics
of the structure (e.g., Bally, 2008). Some examples follow. One recent SN
remnant ($\sim$$10^{5}$ yrs) was identified by Ochsendorf et al. (2015) and
considered at the origin of the H$\alpha$-structure known as Barnard’s loop
(see right panel in Fig. 1). Foley et al. (2023) also considered a newly
discovered cavity in the 3D dust map from Leike et al. (2020), called the
Orion shell (see Fig. 2), which was probably generated by a couple of SN
explosions in the last few Myrs.

Figure 5:  RGB composite image showing: the 3D dust map integrated between 300
and 500 pc in blue colours and contours as in Fig. 2; the 26Al emission map
from Diehl (2002) in red colours and contours. Red contours are between
$2.8\times 10^{-6}$ and $2.4\times 10^{-4}$ cm-2 sr-1 s-1 in steps of
$2.7\times 10^{-5}$ cm-2 sr-1 s-1. The Orion-Taurus ridge is shown with white
contours as in Fig. 1; the $N_{\rm H_{2}}$ map from Kalberla et al. (2020) is
in green colours and contours between $6.3\times 10^{20}$ and $10^{22}$ cm-2
in $\log_{10}$-steps of 0.2. The yellow circles are the same as in Fig. 2.

This hypothesis is supported by the presence of strong 26Al emission in the
region, which has a half-life of 0.7 Myr mostly traces recent feedback events
(e.g., Foley et al., 2023). In Fig. 5, we show the 3D dust map integrated
between 300 and 500 pc combined with the 26Al map based on COMPTEL
$\gamma$-ray data at 1.809 MeV (Diehl, 2002) also presented in Foley et al.
(2023). The image clearly reveals the 26Al emission (in red) within the cavity
seen in dust extinction, encompassed by both the Orion-Taurus ridge (white
contours) and the molecular gas traced by the $N_{\rm H_{2}}$ map (in green)
from Kalberla et al. (2020) smoothed at the resolution of $w_{0}$.

Compression by multiple SN episodes may be at the origin of the organized
structure of $\@vec{B}_{\perp}$ in the arc and of its large strength. Similar
values of $B_{\perp}$ in the superbubble (of the order of tens of $\mu$G) were
also found along its southern Galactic edge (Soler et al., 2018; Joubaud et
al., 2019). Moreover, as shown in Appendix B, using a $B_{\perp}$ value of a
few tens of $\mu$G ($>30-40\,\mu$G) and CRe energy spectra that account for
the $\beta$ distribution in the Orion-Taurus ridge, we were able to recover
the same order of magnitude of $T_{\rm b}$(74 MHz) and $PI$(30 GHz) as
measured by LWA and Planck. Their median values and corresponding errors are
$2534^{+166}_{-210}$ K and $(1.5\pm 1)\times 10^{-5}$ MJy sr-1, respectively
(see Figs. 1 and 10).

##  5 Summary and discussion

In this work, we presented the first detection of diffuse synchrotron emission
from cold-neutral, dusty, and partly molecular gas. We combined data from the
LWA telescope (Dowell et al., 2017), the Planck satellite (Planck 2018
results. I, 2020), and the 3D dust maps of Lallement et al. (2019).

We discovered the physical correlation between one radio loop – the Orion-
Taurus ridge – seen both in total intensity below 100 MHz and in polarized
intensity at 30 GHz, with the density structure of the northern edge of the
OES as traced by interstellar dust extinction. The multi-wavelength
correlation extends over several tens of degrees on the sky ($\sim
25^{\circ}$) with dust located between 300 and 500 pc from the Sun.
Considering an average distance of 400 pc, the Orion-Taurus ridge is about
$\sim 180$ pc wide, which is compatible with the expected edge-size of
superbubbles from SN feedback (e.g., Ntormousi et al., 2017; Grudić et al.,
2021). We note that, in addition to being consistent with the dusty edge of
the OES (Ochsendorf et al., 2015; Pon et al., 2016), the Orion-Taurus ridge
could also be associated with the recently-discovered Orion shell (Foley et
al., 2023).

Whether the Orion-Taurus ridge is linked to the OES or to the Orion shell, it
is most likely still the result of synchrotron radiation from CRe in
relatively dense shock-compressed gas ($N_{\rm H_{2}}\gtrsim 10^{21}$ cm-2)
produced by several episodes of recent SN remnants. This scenario is supported
by enhanced 26Al emission seen in the proximity of the Orion-Taurus ridge
(Diehl, 2002; Foley et al., 2023) and by our analysis of the spectral index of
synchrotron emissivity – $\beta$ – between 50 and 80 MHz from the LWA. After
correcting for possible contamination of free-free absorption along the
sightline, we found a rather flat distribution of $\beta$ in the Orion-Taurus
ridge with a median value equal to $-2.24^{+0.03}_{-0.02}$. Using the
parametric model of Padovani et al. (2018), we explained such flat $\beta$
values by damping the low-energy part of the CRe spectrum below 1 GeV. This
phenomenological fact may be associated with the expected drop in acceleration
efficiency of CRe below a few GeV following the adiabatic phase of SN remnants
(Morlino & Celli, 2021). Because of strong magnetic-field amplification caused
by compression (Ferriere et al., 1991), SN feedback could also explain the
organized structure of the plane-of-the-sky magnetic field (see Fig. 1) and
its strength, $B_{\perp}$, that we derived from the spectral-index analysis of
the order of a few tens of $\mu$G ($>30-40\,\mu$G). These values are
consistent with those already estimated in the southern edge of the OES (e.g.,
Soler et al., 2018) and with the amount of total and polarized synchrotron
emissions from the LWA telescope and the Planck satellite, respectively.

This work highlights the link between the magnetized and cold ISM by means of
low-frequency diffuse synchrotron emission, extending the studies that link
radio-interferometric data in polarization below 200 MHz with tracers of the
neutral ISM, such as HI spectroscopic data, interstellar dust extinction, and
dust thermal polarized emission (Van Eck et al., 2017; Jelić et al., 2018;
Bracco et al., 2020; Turić et al., 2021; Bracco et al., 2022). However, the
case of the Orion-Taurus ridge is the first involving synchrotron total
intensity. Possibly, this is due to the vicinity of the OES to the solar
neighborhood, to the ongoing high-mass star formation in the region, and to
its Galactic location with very low background emission. The advent of future
high-sensitivity facilities in the radio spectrum, such as the C-Band All-Sky
Survey (Jones et al., 2018) and the Square Kilometre Array (Dewdney et al.,
2009), will likely begin revealing structures such as the Orion-Taurus ridge
in a more systematic way. This will open a new window in our understanding of
the complex relationship between cosmic rays, magnetic fields, and
interstellar gas over a wide range of Galactic environments (Dickinson et al.,
2015; Padovani & Galli, 2018).

###### Acknowledgements.

We thank the anonymous referee for the constructive report that improved this
manuscript. We are grateful to F. Mertens for introducing us to the LWA data.
We warmly thank D. Galli, G. Morlino, D. Arzoumanian, I. Grenier, and F.
Boulanger for insightful discussions and comments. We appreciated the kindness
of M. Foley, J. Alves, and R. Diehl for sharing with us the COMPTEL 26Al map
(Diehl, 2002). This work would not have been possible without the inspiring
babbling of Elio and the tireless patience of Delphine. A.B. acknowledges
support from the European Research Council through the Advanced Grant MIST
(FP7/2017-2022, No. 742719). Some of the results in this paper have been
derived using the healpy and HEALPix packages. JDS acknowledges funding by the
European Research Council via the ERC Synergy Grant “ECOGAL – Understanding
our Galactic ecosystem: From the disk of the Milky Way to the formation sites
of stars and planets” (project ID 855130). In the analysis we made use of
astropy (Astropy Collaboration et al., 2018), scipy (Virtanen et al., 2020),
and numpy (Harris et al., 2020).

## References

  * Ackermann et al. (2010) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2010, Phys. Rev. D, 82, 092004
  * Adriani et al. (2011) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, ApJ, 742, 102
  * Aguilar et al. (2014) Aguilar, M., Aisa, D., Alvino, A., et al. 2014, Phys. Rev. Lett., 113, 121102
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bally (2008) Bally, J. 2008, in Handbook of Star Forming Regions, Volume I, ed. B. Reipurth, Vol. 4, 459
  * Berkhuijsen (1971) Berkhuijsen, E. M. 1971, A&A, 14, 359
  * Berkhuijsen et al. (1971) Berkhuijsen, E. M., Haslam, C. G. T., & Salter, C. J. 1971, A&A, 14, 252
  * Borka et al. (2008) Borka, V., Milogradov-Turin, J., & Urošević, D. 2008, Astronomische Nachrichten, 329, 397
  * Bracco et al. (2020) Bracco, A., Jelić, V., Marchal, A., et al. 2020, A&A, 644, L3
  * Bracco et al. (2022) Bracco, A., Ntormousi, E., Jelić, V., et al. 2022, A&A, 663, A37
  * Cabral & Leedom (1993) Cabral, B. & Leedom, L. C. 1993, in Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’93 (New York, NY, USA: Association for Computing Machinery), 263–270
  * Cummings et al. (2016) Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18
  * Dewdney et al. (2009) Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proceedings, 97, 1482
  * Dickinson et al. (2015) Dickinson, C., Beck, R., Crocker, R., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 102
  * Diehl (2002) Diehl, R. 2002, New A Rev., 46, 547
  * Dowell et al. (2017) Dowell, J., Taylor, G. B., Schinzel, F. K., Kassim, N. E., & Stovall, K. 2017, MNRAS, 469, 4537
  * Ferrière (2020) Ferrière, K. 2020, Plasma Physics and Controlled Fusion, 62, 014014
  * Ferriere et al. (1991) Ferriere, K. M., Mac Low, M.-M., & Zweibel, E. G. 1991, ApJ, 375, 239
  * Ferriere & Zweibel (1991) Ferriere, K. M. & Zweibel, E. G. 1991, ApJ, 383, 602
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, ApJS, 146, 407
  * Foley et al. (2023) Foley, M. M., Goodman, A., Zucker, C., et al. 2023, ApJ, 947, 66
  * Gaensler et al. (2011) Gaensler, B. M., Haverkorn, M., Burkhart, B., et al. 2011, Nature, 478, 214
  * Ginzburg & Syrovatskii (1965) Ginzburg, V. L. & Syrovatskii, S. I. 1965, ARA&A, 3, 297
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Großschedl et al. (2021) Großschedl, J. E., Alves, J., Meingast, S., & Herbst-Kiss, G. 2021, A&A, 647, A91
  * Grudić et al. (2021) Grudić, M. Y., Kruijssen, J. M. D., Faucher-Giguère, C.-A., et al. 2021, MNRAS, 506, 3239
  * Guzmán et al. (2011) Guzmán, A. E., May, J., Alvarez, H., & Maeda, K. 2011, A&A, 525, A138
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Haslam et al. (1971) Haslam, C. G. T., Kahn, F. D., & Meaburn, J. 1971, A&A, 12, 388
  * Heiles (1976) Heiles, C. 1976, ApJ, 208, L137
  * Heiles (1997) Heiles, C. 1997, ApJS, 111, 245
  * Jelić et al. (2018) Jelić, V., Prelogović, D., Haverkorn, M., Remeijn, J., & Klindžić, D. 2018, A&A, 615, L3
  * Jones et al. (2018) Jones, M. E., Taylor, A. C., Aich, M., et al. 2018, MNRAS, 480, 3224
  * Joubaud et al. (2019) Joubaud, T., Grenier, I. A., Ballet, J., & Soler, J. D. 2019, A&A, 631, A52
  * Joubaud et al. (2020) Joubaud, T., Grenier, I. A., Casandjian, J. M., Tolksdorf, T., & Schlickeiser, R. 2020, A&A, 635, A96
  * Kalberla et al. (2020) Kalberla, P. M. W., Kerp, J., & Haud, U. 2020, A&A, 639, A26
  * Kassim (1989) Kassim, N. E. 1989, ApJ, 347, 915
  * Kataoka et al. (2021) Kataoka, J., Yamamoto, M., Nakamura, Y., et al. 2021, ApJ, 908, 14
  * Lallement et al. (2019) Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135
  * Leike et al. (2020) Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138
  * Martire et al. (2023) Martire, F. A., Banday, A. J., Martínez-González, E., & Barreiro, R. B. 2023, J. Cosmology Astropart. Phys., 2023, 049
  * Mertsch & Sarkar (2013) Mertsch, P. & Sarkar, S. 2013, J. Cosmology Astropart. Phys., 2013, 041
  * Milogradov-Turin & Urošević (1997) Milogradov-Turin, J. & Urošević, D. 1997, Bulletin Astronomique de Belgrade, 155, 41
  * Miville-Deschênes et al. (2008) Miville-Deschênes, M. A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093
  * Morlino & Celli (2021) Morlino, G. & Celli, S. 2021, MNRAS, 508, 6142
  * Ntormousi et al. (2017) Ntormousi, E., Dawson, J. R., Hennebelle, P., & Fierlinger, K. 2017, A&A, 599, A94
  * Ochsendorf et al. (2015) Ochsendorf, B. B., Brown, A. G. A., Bally, J., & Tielens, A. G. G. M. 2015, ApJ, 808, 111
  * Orlando (2018) Orlando, E. 2018, MNRAS, 475, 2724
  * Padovani et al. (2021) Padovani, M., Bracco, A., Jelić, V., Galli, D., & Bellomi, E. 2021, A&A, 651, A116
  * Padovani & Galli (2018) Padovani, M. & Galli, D. 2018, A&A, 620, L4
  * Padovani et al. (2018) Padovani, M., Galli, D., Ivlev, A. V., Caselli, P., & Ferrara, A. 2018, A&A, 619, A144
  * Panopoulou et al. (2021) Panopoulou, G. V., Dickinson, C., Readhead, A. C. S., Pearson, T. J., & Peel, M. W. 2021, ApJ, 922, 210
  * Planck 2013 results. XI. (2014) Planck 2013 results. XI. 2014, A&A, 571, A11
  * Planck 2015 results. X. (2016) Planck 2015 results. X. 2016, A&A, 594, A10
  * Planck 2015 results. XXV. (2016) Planck 2015 results. XXV. 2016, A&A, 594, A25
  * Planck 2018 results. I (2020) Planck 2018 results. I. 2020, A&A, 641, A1
  * Planck Collaboration Int. XXXV (2016) Planck Collaboration Int. XXXV. 2016, A&A, 586, A138
  * Pon et al. (2016) Pon, A., Ochsendorf, B. B., Alves, J., et al. 2016, ApJ, 827, 42
  * Reynolds & Ogden (1979) Reynolds, R. J. & Ogden, P. M. 1979, ApJ, 229, 942
  * Reynolds et al. (1998) Reynolds, R. J., Tufte, S. L., Haffner, L. M., Jaehnig, K., & Percival, J. W. 1998, PASA, 15, 14
  * Rybicki & Lightman (1979) Rybicki, G. B. & Lightman, A. P. 1979, Radiative processes in astrophysics
  * Salter (1983) Salter, C. J. 1983, Bulletin of the Astronomical Society of India, 11, 1
  * Sofue (2015) Sofue, Y. 2015, MNRAS, 447, 3824
  * Soler et al. (2018) Soler, J. D., Bracco, A., & Pon, A. 2018, A&A, 609, L3
  * Stanislavsky et al. (2023) Stanislavsky, L. A., Bubnov, I. N., Konovalenko, A. A., Stanislavsky, A. A., & Yerin, S. N. 2023, A&A, 670, A157
  * Stone et al. (2019) Stone, E. C., Cummings, A. C., Heikkila, B. C., & Lal, N. 2019, Nature Astronomy, 3, 1013
  * Tahani et al. (2022) Tahani, M., Glover, J., Lupypciw, W., et al. 2022, A&A, 660, L7
  * Taylor et al. (2012) Taylor, G. B., Ellingson, S. W., Kassim, N. E., et al. 2012, Journal of Astronomical Instrumentation, 1, 1250004
  * Thomson et al. (2021) Thomson, A. J. M., Landecker, T. L., McClure-Griffiths, N. M., et al. 2021, MNRAS, 507, 3495
  * Turić et al. (2021) Turić, L., Jelić, V., Jaspers, R., et al. 2021, A&A, 654, A5
  * Van Eck et al. (2017) Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2017, A&A, 597, A98
  * Vidal et al. (2015) Vidal, M., Dickinson, C., Davies, R. D., & Leahy, J. P. 2015, MNRAS, 452, 656
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298

##  Appendix A Distinct structures in the 3D dust map

The Orion-Taurus ridge under study only shows morphological correlation with
the 3D dust map at the position of the northern edge of the OES ($\approx 400$
pc).

Figure 6: Same as in Fig. 2 but for two different distance ranges, namely
between 100 and 300 pc (top panel) and between 500 and 700 pc (bottom panel).

Figure 7: Same as in the left panel of Fig. 3 but for two distinct distance
ranges.

In this section we show the visual comparison (Fig. 6) and the $r^{\rm MAX}$
distributions (Fig. 7) between the Orion-Taurus ridge defined at 74 MHz and
the 3D dust map integrated at different locations along the LOS: between 100
and 300 pc and between 500 and 700 pc. The structures in the maps are
uncorrelated with $\Delta r^{\rm MAX}_{\rm m}=14^{+2}_{-17}$ ∘ and $\Delta
r^{\rm MAX}_{\rm m}=-3^{+14}_{-8}$ ∘, respectively.

In the bottom panel of Fig. 7 the dust ring reported in Tahani et al. (2022)
can be seen in the proximity of the OMC.

##  Appendix B Cosmic-ray electron energy spectrum

In this Appendix we summarize and detail the parametric model for the energy
flux of CRe – $j_{e}(E)$ – from Padovani et al. (2018) that was used in this
work to propose a possible mechanism for the $\beta$ distribution of
synchrotron emissivity in the Orion-Taurus ridge (see Sect. 4).

###  B.1 Explaining the $\beta$ distribution

The synchrotron $T_{\rm b}(\nu)$ is proportional to the integral along the LOS
of the sum of the specific emissivities polarized along and across
$\@vec{B_{\perp}}$, namely $\epsilon_{\nu,\parallel}$ and
$\epsilon_{\nu,\perp}$ (Ginzburg & Syrovatskii 1965; Padovani et al. 2021;
Bracco et al. 2022). These two quantities depend upon the CRe energy spectrum,
given by $j_{e}(E)$, and the value of $B_{\perp}(\@vec{r})$ at position
$\@vec{r}$. Changing $j_{e}(E)$ and $B_{\perp}(\@vec{r})$, one can control the
frequency dependence of $T_{\rm b}(\nu)$, thus the value of $\beta$.

Figure 8: Parametric models for the CRe energy spectrum from Padovani et al.
(2018). See Eq. 2 for reference.

Using state-of-the-art magneto-hydrodynamical numerical simulations of the
multiphase ISM, Padovani et al. (2021) thoroughly explored the properties of
$\beta$ as a function of $B_{\perp}(\@vec{r})$ and $j_{e}(E)$. They concluded
that the shape of $j_{e}(E)$ has a particularly strong impact on actual
observables of synchrotron emissivity at low frequency ($\lesssim$ 10 GHz). In
their work the authors mostly relied on the local $j_{e}(E)$ deduced in
Padovani et al. (2018) from real data of interstellar CRe with $E$ between a
few hundred MeV (Voyager 1-2, Cummings et al. 2016; Stone et al. 2019) and
hundreds of GeV (Fermi-LAT, Pamela, and AMS, Ackermann et al. 2010; Adriani et
al. 2011; Aguilar et al. 2014, respectively).

Figure 9: 2D correlation plots between $\beta$ and $a$ as a function of the
averaged $B_{\perp}$ along the sight line, $\langle B_{\perp}\rangle$ (in
colours and contours). In the top panel $E_{0}=710$ MeV; in the bottom panel,
$E_{0}=1$ GeV. The grey-shaded areas and the dahed-black lines are $\beta^{\rm
c}_{\rm m}=-2.24^{+0.03}_{-0.02}$ (see Sect. 4).

We noticed that by applying the same $j_{e}(E)$ as in Padovani et al. (2018) –
hereafter, $j_{e0}$ – we could not reproduce the $\beta$ distribution of the
Orion-Taurus ridge measured below 100 MHz with the LWA. In Padovani et al.
(2018), $j_{e0}$ is parameterized as follows:

$j_{e0}(E)=C\frac{E^{a}}{(E+E_{0})^{(3.2+a)}},$ (2)

where $C=2.1\times 10^{18}$ eV-1 s-1 cm-2 sr-1, $a$ $=$ $-1.3$, and $E_{0}$
$=$ $710$ MeV.

In this work, we explored the possibility of changing $a$ and $E_{0}$ in order
to obtain values of $\beta$ consistent with $\beta^{\rm c}_{\rm m}$ (see Sect.
4). In Fig. 8, we illustrate the $j_{e}(E)$ variation as a function of $a$ and
$E_{0}$ compared to $j_{e0}$ (see black circles). While $a$ impacts the slope
of $j_{e}(E)$ at low energy, $E_{0}$ may also slightly affect the high-energy
end of $j_{e}(E)$. Using such analytical models, as shown in Fig. 9, we
explored a wide range of parameters. We could reproduce the flat $\beta$
observed between 50 and 80 MHz (see the grey horizontal band in Fig. 9) by
only considering values of $a\gtrsim-0.5$ (lower population of CRe below a few
hundred MeV compared to $j_{e0}$) and averaged $B_{\perp}$ along the sight
line, $\langle B_{\perp}\rangle$, larger than 30 $\mu$G. This conclusion holds
for both $E_{0}=710$ MeV and $E_{0}=1$ GeV. We note that we could also reach
$\beta$ values consistent with $\beta^{\rm c}_{\rm m}$ even for smaller values
of $a$ and $\langle B_{\perp}\rangle$ if only $E_{0}$ was larger. Looking at
Fig. 8, this would correspond to having $j_{e}$ flatter than $j_{e0}$ at high
energy too. As explained by Morlino & Celli (2021), this effect could also be
expected in SN remnants at high energies larger than 10 TeV, where a
flattening of the CRe energy spectrum is predicted as the result of electrons
escaping from the shock precursors. However, high-energy Fermi-LAT
observations suggest that the CR-proton energy spectrum of the OES does not
show significant deviations from the local interstellar one (i.e., Joubaud et
al. 2020). If this constraint holds for CRe at high energy, the CRe energy
spectrum at low energies ($<$ GeV) may be different due to, for example, large
energy losses caused by the enhanced magnetic field in OES (Morlino & Celli
2021). A more detailed analysis of the high-energy CRe is left to follow-up
work since their impact on the low-frequency radio emission is rather
negligible (e.g. Fig. 3 in Padovani et al. 2021).

Figure 10: Estimates of $T_{\rm b}$(74 MHz) – top – and $PI_{\rm MAX}$(30 GHz)
– bottom – of the Orion-Taurus ridge from models of Padovani et al. (2018,
2021). The quantities are derived as a function of $\langle B_{\perp}\rangle$,
$a$, and $E_{0}$. The horizontal black dashed lines, with grey errors,
represent the corresponding median values reported in Sect. 4.

###  B.2 Synchrotron intensity estimates

Our results allowed us to give order-of-magnitude estimates of the synchrotron
emission from the Orion-Taurus ridge that we could compare with LWA and Planck
observations (see Sec. 2). We refer to the analytical formulation described by
Eqs. (2-7) of Padovani et al. (2021), which explains in detail how $T_{\rm b}$
and $PI$ can be derived from $\epsilon_{\nu,\perp}$ and
$\epsilon_{\nu,\parallel}$ given $B_{\perp}$, $j_{e}$, and the observational
frequency $\nu$ (see also, Ginzburg & Syrovatskii 1965). In short, assuming a
homogeneous distribution of CRe and magnetic fields, we made use of the
following relations, $T_{\rm b}$ $\propto$
$\langle\epsilon_{\nu,\parallel}+\epsilon_{\nu,\perp}\rangle$ and
$PI\propto\langle\epsilon_{\nu,\perp}-\epsilon_{\nu,\parallel}\rangle$, where
brackets indicate averages along the LOS. For consistency with the 3D dust
map, we considered a LOS path-length ($L$) of 200 pc. The estimates that we
obtained represent upper limits to $T_{\rm b}$ and $PI$ at 74 MHz and 30 GHz,
respectively, as most likely the Orion-Taurus ridge is not uniformly
distributed along $L$.

We also note that in the case of $PI$ we provided maximum polarization
estimates ($PI_{\rm MAX}$) without including realistic LOS depolarization
effects (Miville-Deschênes et al. 2008; Gaensler et al. 2011). The tangling of
the magnetic field lines along the LOS and within the telescope beam reduces
the amount of observed polarization. Depending on the LOS depth and on the
turbulent component of the magnetic field, the drop in the observed
polarization can range from a fraction to tens of percent for a LOS extent of
$L$ $=$ $200$ pc (see Fig. 6 in Miville-Deschênes et al. 2008). In our case,
however, the polarization fraction ($PI/I$, where $I$ is the total intensity)
of all synchrotron models was the intrinsic theoretical value of 75% (Rybicki
& Lightman 1979). We neglected the effect of Faraday rotation, which is
unimportant at 30 GHz (Ferrière 2020).

Figure 10 shows our estimates of $T_{\rm b}$(74 MHz) and $PI_{\rm MAX}$(30
GHz) after applying a beam of $2^{\circ}$ FWHM resolution. Despite all caveats
listed above, we have shown that given magnetic-field strengths of a few tens
of $\mu$G ($>40\,\mu$G) and the CRe energy properties of the Orion-Taurus
ridge derived in Appendix B.1, we could retrieve the right orders of magnitude
of both $T_{\rm b}$(74 MHz) and $PI$(30 GHz) from LWA and Planck,
respectively.