Magnetic turbulence spectra and intermittency in the heliosheath and in the local interstellar medium

TL;DR

This study analyzes magnetic turbulence spectra and intermittency in the heliosheath and local interstellar medium using Voyager data, revealing distinct turbulence regimes, anisotropy, and intermittency characteristics across scales.

Contribution

It provides the first highly resolved statistical analysis of magnetic turbulence across a broad range of scales in the heliosheath and LISM, including the MHD-to-kinetic transition.

Findings

Identification of energy-injection and inertial-cascade regimes in the heliosheath.

Observation of anisotropic and compressive turbulence modes.

Intermittency occurs only at small scales in the interstellar medium.

Abstract

The understanding of inertial-scale dynamics in the heliosheath is not yet thorough. Magnetic field fluctuations across the inner heliosheath (IHS) and the local interstellar medium (LISM) are here considered to provide accurate and highly resolved statistics over different plasma conditions between 88 and 136 au. By using the unique in situ 48 s measurements from the Voyager Interstellar Mission, we investigate different fluctuation regimes at the magnetohydrodynamic (MHD) scales, down to the MHD-to-kinetic transition. We focus on a range of scales exceeding five frequency decades (5 x 10-8 < f < 10-2 Hz), which is unprecedented in literature analysis. A set of magnetic field data for eight intervals in the IHS, in both unipolar and sector regions, and four intervals in the LISM is used for the analysis. Results are set forth in terms of the power spectral density, spectral compressibility, structure functions, and intermittency of magnetic field increments. In the heliosheath, we identify the energy-injection regime displaying a ~1/f energy decay, and the inertial-cascade regime. Here, the power spectrum is anisotropic and dominated by compressive modes, with intermittency that can reach kurtosis values of up to 10. In the interstellar medium the structure of turbulence is anisotropic as well, with transverse fluctuations clearly prevailing after 2015 May. Here, we show that intermittent features occur only at scales smaller than 10-6 Hz.

Figures

Tables

Table 1. Table 1: Averaged quantities at Voyager 2 in the IHS. Plasma quantities are computed from PLS data available in the NASA COHO web site. The table reports r SC subscript 𝑟 SC r_{\mathrm{SC}} , the Sun-V1 distance; V 0 = ( V 0 R 2 + V 0 T 2 + V 0 N 2 ) 1 / 2 subscript 𝑉 0 superscript superscript subscript 𝑉 subscript 0 𝑅 2 superscript subscript 𝑉 subscript 0 𝑇 2 superscript subscript 𝑉 subscript 0 𝑁 2 1 2 V_{0}=(V_{0_{R}}^{2}+V_{0_{T}}^{2}+V_{0_{N}}^{2})^{1/2} and B 0 subscript 𝐵 0 B_{0} , the magnitudes of the average velocity and magnetic field, respectively; B 𝐵 B , the average magnetic field strength. For the thermal protons: f cp subscript 𝑓 cp f_{\mathrm{cp}} , the cyclotron frequency; n p subscript 𝑛 p n_{\mathrm{p}} , the average density; T p subscript 𝑇 p T_{\mathrm{p}} , average temperature; V A subscript 𝑉 A V_{\mathrm{A}} , the Alfvén velocity ; β p = ⟨ 2 ​ μ 0 ​ n p ​ k B ​ T p / B 2 ⟩ subscript 𝛽 p delimited-⟨⟩ 2 subscript 𝜇 0 subscript 𝑛 p subscript 𝑘 B subscript 𝑇 p superscript 𝐵 2 \beta_{\mathrm{p}}=\langle 2\mu_{0}n_{\mathrm{p}}k_{\mathrm{B}}T_{\mathrm{p}}/B^{2}\rangle , the beta of thermal ions; β p ​ 1 ​ keV subscript 𝛽 p 1 keV \beta_{\mathrm{p\ 1keV}} , the beta of 1-keV PUIs ( n 1 ​ k ​ e ​ V ≈ 0.2 ​ n subscript 𝑛 1 𝑘 𝑒 𝑉 0.2 𝑛 n_{1keV}\approx 0.2~{}n ); r ip subscript 𝑟 ip r_{\mathrm{ip}} , the ion inertial radius; r cp subscript 𝑟 cp r_{\mathrm{cp}} , the gyroradius. r cp ​ 1 ​ keV subscript 𝑟 cp 1 keV r_{\mathrm{cp\ 1keV}} is the ion gyroradius of a 1-keV pickup proton. Via the Taylor’s approximation, frequencies are also converted in the spacecraft reference frame f SC ≈ V SW / ( 2 ​ r ) subscript 𝑓 SC subscript 𝑉 SW 2 𝑟 f_{\mathrm{SC}}\approx V_{\mathrm{SW}}/(2r) . Eventually, the 1-eddy-turnover frequency f e subscript 𝑓 𝑒 f_{e} is shown (Eq. 7 ).

Voyager 2	SHS1	UHS1	UHS2	SHS2
$r_{\mathrm{SC}}\left(\mathrm{AU}\right)$	89.0	94.1	100.2	104.9
$V_0\left(\mathrm{km}{s}^{-1}\right)$	157	151	154	153
$B_0\left(\mathrm{nT}\right)$	0.062	0.072	0.090	0.030
$B\left(\mathrm{nT}\right)$	0.096	0.086	0.128	0.100
$n_p\left({\mathrm{cm}}^{-3}\right)$	1.1$\times {10}^{-3}$	1.0$\times {10}^{-3}$	1.9$\times {10}^{-3}$	1.8$\times {10}^{-3}$
$T_p\left({10}^{4}K\right)$	6.35	4.47	5.63	5.17
$V_A\left(\mathrm{km}{s}^{-1}\right)$	63.4	59.2	63.9	51.6
${\beta }_p$	0.54	0.46	0.54	1.00
${\beta }_{p1keV}$	19.6	24.3	22.2	45.0
$r_{\mathrm{ip}}\left(\mathrm{km}\right)$	6852	7186	5214	5357
$r_{\mathrm{cp}}\left(\mathrm{km}\right)$	3502	3301	2485	3038
$r_{\mathrm{cp}1\mathrm{keV}}\left(\mathrm{km}\right)$	53000	53165	35695	45480
$f_{\mathrm{ip}\mathrm{SC}}\left(\mathrm{mHz}\right)$	11.5	10.5	14.8	14.3
$f_{\mathrm{cp}}\left(\mathrm{mHz}\right)$	1.47	1.31	1.95	1.15
$f_{\mathrm{cp}\mathrm{SC}}\left(\mathrm{mHz}\right)$	22.4	22.8	30.9	25.2
$f_{\mathrm{cp}1\mathrm{keV}\mathrm{SC}}\left(\mathrm{mHz}\right)$	1.65	1.42	2.16	1.68
$f_e\left(\mathrm{Hz}\right)$	1.6$\times {10}^{-6}$	8$\times {10}^{-7}$	4$\times {10}^{-7}$	4.5$\times {10}^{-7}$

Table 2. Table 2: Magnetic field fluctuations properties at V2 in the IHS. All quantities are defined in § 3.1 . Spectral breaks and indexes of perpendicular fluctuations refer to the total power spectrum P ​ [ B ⟂ ] = P ​ [ B ⟂ 1 ] + P ​ [ B ⟂ 2 ] P delimited-[] subscript 𝐵 perpendicular-to P delimited-[] subscript 𝐵 perpendicular-to absent 1 P delimited-[] subscript 𝐵 perpendicular-to absent 2 \mathrm{P}[B_{\perp}]=\mathrm{P}[B_{\perp 1}]+\mathrm{P}[B_{\perp 2}] .

Voyager 2	SHS1	UHS1	UHS2	SHS2
$E_m\left({\mathrm{nT}}^2\right)$	6.76$\times {10}^{-3}$	3.42$\times {10}^{-3}$	1.17$\times {10}^{-2}$	1.18$\times {10}^{-2}$
$C_2$	0.29	0.42	0.35	0.43
$C_2$	0.35	0.63	0.48	0.64
$I_{\parallel }$	0.50	0.45	0.52	2.29
$I_{⟂1}$	0.62	0.37	0.50	2.14
$I_{⟂2}$	0.44	0.25	0.49	1.06
$I$	1.04	0.74	1.00	3.70
$\delta B_{\mathrm{mv}}$	0.051	0.043	0.067	0.088
${\theta }_{mv}$	106∘	153∘	86∘	51∘
$f_{\mathrm{b1}}$ (Hz)	5$\times {10}^{-6}$	${10}^{-5}$	$7\times {10}^{-6}$	$7\times {10}^{-7}$
$f_{\mathrm{b2}}$ (Hz)	2$\times {10}^{-4}$	$\cdots$	${10}^{-4}$	7$\times {10}^{-5}$
${\alpha }_1$	-0.98$\pm$0.10	-1.10$\pm$0.03	-1.18$\pm$0.05	-1.25$\pm$0.13
${\alpha }_2$	-1.64$\pm$0.01	-1.78$\pm$0.06	-1.78$\pm$0.13	-1.58$\pm$0.02
$f_{\mathrm{b1},\parallel }$ (Hz)	2$\times {10}^{-6}$	3$\times {10}^{-6}$	2$\times {10}^{-6}$	7$\times {10}^{-7}$
$f_{\mathrm{b2},\parallel }$ (Hz)	2$\times {10}^{-4}$	${10}^{-4}$	5$\times {10}^{-5}$	6$\times {10}^{-5}$
${\alpha }_{1,\parallel }$	-0.74$\pm$0.10	-0.97$\pm$0.12	-0.84$\pm$0.05	-1.15$\pm$0.02
${\alpha }_{2,\parallel }$	-1.62$\pm$0.02	-1.60$\pm$0.10	-1.34$\pm$0.05	-1.57$\pm$0.03
${\alpha }_{3,\parallel }$	-1.97$\pm$0.02	…	-2.08$\pm$0.02	-2.02$\pm$0.08
$f_{\mathrm{b1},⟂}$ (Hz)	5$\times {10}^{-6}$	2$\times {10}^{-5}$	8$\times {10}^{-6}$	7$\times {10}^{-7}$
$f_{\mathrm{b2},⟂}$ (Hz)	2$\times {10}^{-4}$	$\cdots$	$\cdots$	7$\times {10}^{-5}$
${\alpha }_{1,⟂}$	-0.80$\pm$0.10	-1.25$\pm$0.04	-1.13$\pm$0.03	-1.11$\pm$0.10
${\alpha }_{2,⟂}$	-1.64$\pm$0.02	-1.87$\pm$0.06	-2.05$\pm$ 0.05	-1.58$\pm$0.02
${\alpha }_{3,⟂}$	-1.70$\pm$0.05	$\cdots$	$\cdots$	-1.70$\pm$0.05
${\zeta }_{1,\parallel }$ $\left({\zeta }_{1,\parallel }^{\mathrm{ess}}\right)$	0.33 (0.35)	0.30 (0.33)	0.28 (0.41)	0.40 (0.41)
${\zeta }_{2,\parallel }$ $\left({\zeta }_{2,\parallel }^{\mathrm{ess}}\right)$	0.67 (0.68)	0.57 (0.67)	0.48 (0.74)	0.75 (0.75)
${\zeta }_{3,\parallel }$ $\left({\zeta }_{3,\parallel }^{\mathrm{ess}}\right)$	1.02 (1)	0.82 (1)	0.67 (1)	1.05 (1)
${\zeta }_{4,\parallel }$ $\left({\zeta }_{4,\parallel }^{\mathrm{ess}}\right)$	1.38 (1.30)	1.04 (1.31)	0.89 (1.20)	1.31 (1.18)
${\zeta }_{1,⟂}$ $\left({\zeta }_{1,⟂}^{\mathrm{ess}}\right)$	0.40 (0.37)	0.25 (0.34)	0.37 (0.38)	0.33 (0.36)
${\zeta }_{2,⟂}$ $\left({\zeta }_{2,⟂}^{\mathrm{ess}}\right)$	0.82 (0.71)	0.45 (0.67)	0.62 (0.71)	0.62 (0.70)
${\zeta }_{3,⟂}$ $\left({\zeta }_{3,⟂}^{\mathrm{ess}}\right)$	1.25 (1)	0.64 (1)	0.81 (1)	0.87 (1)
${\zeta }_{4,⟂}$ $\left({\zeta }_{4,⟂}^{\mathrm{ess}}\right)$	1.69 (1.24)	0.82 (1.30)	0.98 (1.26)	1.10 (1.25)

Table 3. Table 3: Averaged quantities at V1 in the IHS. Since the PLS subsystem is not operative, the velocity is derived from LECP and CRS subsystems. Here, we report data from Figure 1 in Krimigis et al. ( 2011 ) and Figure 1 in Richardson & the Voyager Team ( 2016 ) .

Voyager 1	A1	B1	C1	D1
$r_{\mathrm{SC}}\left(\mathrm{AU}\right)$	109.5	115.7	117.9	118.8
$V_0\left(\mathrm{km}{s}^{-1}\right)$	65	40	40	40
$B_0\left(\mathrm{nT}\right)$	0.083	0.132	0.195	0.124
$B\left(\mathrm{nT}\right)$	0.086	0.140	0.203	0.148
$f_{\mathrm{cp}}\left(\mathrm{mHz}\right)$	1.32	2.14	3.10	0.81

Table 4. Table 4: Magnetic field fluctuations properties at V1 in the IHS.

	A1	B1	C1	D1
$E_m\left({\mathrm{nT}}^2\right)$	1.67$\times {10}^{-3}$	5.03$\times {10}^{-3}$	4.50$\times {10}^{-3}$	4.05$\times {10}^{-3}$
$C_1$	0.61	0.55	0.62	0.42
$C_2$	0.95	0.76	0.74	0.62
$I_{\parallel }$	0.34	0.40	0.24	0.28
$I_{⟂1}$	0.11	0.17	0.10	0.27
$I_{⟂2}$	0.13	0.22	0.12	0.13
$I$	0.42	0.54	0.31	0.46
$\delta B_{\mathrm{mv}}$	0.037	0.061	0.069	0.055
${\theta }_{mv}$	166∘	147∘	141∘	42∘
$f_{\mathrm{b1}}$ (Hz)	$2\times {10}^{-5}$	${10}^{-5}$	$2\times {10}^{-5}$	$\cdots$
${\alpha }_1$	-1.17$\pm$0.09	-1.31$\pm$0.15	-1.60$\pm$0.18	$\cdots$
${\alpha }_2$	-2.25$\pm$0.12	-2.31$\pm$0.04	-2.23$\pm$0.09	-1.72$\pm$0.05
$f_{\mathrm{b1}\parallel }$ (Hz)	$2\times {10}^{-5}$	${10}^{-5}$	$2\times {10}^{-5}$	$\cdots$
${\alpha }_{1\parallel }$	-1.18$\pm$0.06	-1.21$\pm$0.07	-1.45$\pm 0.10$	-1.95$\pm 0.10$
${\alpha }_{2\parallel }$	-2.65$\pm$0.10	-2.50$\pm$0.03	-2.52$\pm$0.10	-1.80$\pm 0.05$
$f_{\mathrm{b1}⟂}$ (Hz)	$4\times {10}^{-6}$	$\cdots$	$\cdots$	$\cdots$
${\alpha }_{1⟂}$	-1.35$\pm$0.04	-1.52$\pm$0.05	-1.57$\pm 0.09$	-1.76$\pm 0.10$
${\alpha }_{2⟂}$	-1.59$\pm$0.13	-1.77$\pm$0.06	-1.70$\pm 0.1$	-1.80$\pm 0.05$
${\zeta }_{1,\parallel }$ $\left({\zeta }_{1,\parallel }^{\mathrm{ess}}\right)$	0.48 (0.38)	0.62 (0.36)	0.44 (0.39)	0.36 (0.33)
${\zeta }_{2,\parallel }$ $\left({\zeta }_{2,\parallel }^{\mathrm{ess}}\right)$	0.93 (0.71)	1.19 (0.70)	0.90 (0.72)	0.62 (0.65)
${\zeta }_{3,\parallel }$ $\left({\zeta }_{3,\parallel }^{\mathrm{ess}}\right)$	1.31 (1)	1.63 (1)	1.34 (1)	0.81 (1)
${\zeta }_{4,\parallel }$ $\left({\zeta }_{4,\parallel }^{\mathrm{ess}}\right)$	1.64 (1.24)	1.99 (1.24)	1.72 (1.27)	0.98 (1.23)
${\zeta }_{1,⟂}$ $\left({\zeta }_{1,⟂}^{\mathrm{ess}}\right)$	0.34 (0.36)	0.41 (0.36)	0.43 (0.32)	0.36 (0.38)
${\zeta }_{2,⟂}$ $\left({\zeta }_{2,⟂}^{\mathrm{ess}}\right)$	0.67 (0.70)	0.81 (0.69)	1.02 (0.66)	0.74 (0.70)
${\zeta }_{3,⟂}$ $\left({\zeta }_{3,⟂}^{\mathrm{ess}}\right)$	0.97 (1)	1.17 (1)	0.87 (1)	1.14 (1)
${\zeta }_{4,⟂}$ $\left({\zeta }_{4,⟂}^{\mathrm{ess}}\right)$	1.23 (1.24)	1.48 (1.25)	2.30 (1.35)	1.56 (1.28)

Table 5. Table 5: Magnetic field fluctuations properties in the LISM.

	L1	L2	L3	L4
$E_m\left({\mathrm{nT}}^2\right)$	2.84$\times {10}^{-4}$	2.38$\times {10}^{-4}$	3.72$\times {10}^{-4}$	6.86$\times {10}^{-4}$
$C_1$	0.30	0.40	0.48	0.12
$C_2$	0.41	0.54	0.65	0.12
$I_{\parallel }$	0.014	0.018	0.023	0.015
$I_{⟂1}$	0.019	0.016	0.014	0.030
$I_{⟂2}$	0.009	0.012	0.014	0.038
$I$	0.029	0.030	0.035	0.057
$\delta B_{\mathrm{mv}}$	0.013	0.011	0.015	0.019
${\theta }_{mv}$	65∘	23∘	173∘	98∘
$\alpha$	-1.57$\pm$0.05	-1.65$\pm$0.10	-1.55$\pm$0.10	-1.77$\pm$0.10
${\alpha }_{\parallel }$	-1.59$\pm$0.05	-1.60$\pm$0.10	-1.57$\pm$0.10	-1.40$\pm$0.06
${\alpha }_{⟂}$	-1.54$\pm$0.02	-1.60$\pm$0.05	-1.54$\pm$0.10	-1.90$\pm$0.05
${\zeta }_{1,\parallel }$ $\left({\zeta }_{1,\parallel }^{\mathrm{ess}}\right)$	0.18 (0.33)	0.15 (0.34)	0.14 (0.38)	0.15 (0.35)
${\zeta }_{2,\parallel }$ $\left({\zeta }_{2,\parallel }^{\mathrm{ess}}\right)$	0.36 (0.67)	0.28 (0.68)	0.27 (0.71)	0.29 (0.68)
${\zeta }_{3,\parallel }$ $\left({\zeta }_{3,\parallel }^{\mathrm{ess}}\right)$	0.55 (1)	0.40 (1)	0.39 (1)	0.43 (1)
${\zeta }_{4,\parallel }$ $\left({\zeta }_{4,\parallel }^{\mathrm{ess}}\right)$	0.73 (1.32)	0.51 (1.31)	0.49 (1.25)	0.56 (1.30)
${\zeta }_{1,⟂}$ $\left({\zeta }_{1,⟂}^{\mathrm{ess}}\right)$	0.35 (0.36)	0.25 (0.38)	0.21 (0.40)	0.40 (0.43)
${\zeta }_{2,⟂}$ $\left({\zeta }_{2,⟂}^{\mathrm{ess}}\right)$	0.67 (0.71)	0.48 (0.71)	0.37 (0.72)	0.72 (0.76)
${\zeta }_{3,⟂}$ $\left({\zeta }_{3,⟂}^{\mathrm{ess}}\right)$	0.96 (1)	0.69 (1)	0.5 (1)	0.96 (1)
${\zeta }_{4,⟂}$ $\left({\zeta }_{4,⟂}^{\mathrm{ess}}\right)$	1.21 (1.23)	0.88 (1.24)	0.59 (1.24)	1.12 (1.17)

Equations

$$C_1=⟨{\left(\frac{{\mathbf{B}}_0\cdot \mathbf{\delta }\mathbf{B}}{B_0\delta B}\right)}^2⟩,C_2=\frac{⟨(\mathbf{b}\cdot \mathbf{\delta }\mathbf{B}{)}^2⟩}{⟨\delta B^2⟩}.$$

$$I=⟨\left|\frac{\delta B}{B_0}\right|⟩,I_j=⟨\left|\frac{\delta B_j}{B_0}\right|⟩,j=\{\parallel ,{\perp }_1,{\perp }_2\},$$

$$S_{p,j}(\tau )=⟨|B_j(t)-B_j(t+\tau ){|}^p⟩,j=\{\parallel ,{\perp }_1,{\perp }_2\}.$$

$$S_{p,j}({\tau }_k)=\frac{1}{N({\tau }_k)}\sum _{i=1}^{N({\tau }_k)}|B_j(t_i)-B_j(t_i+{\tau }_k){|}^p,$$

$${\tau }_k=k\cdot \Delta t_{s}k=1,\dots ,n$$

$$K(\tau )=\frac{S_4(\tau )}{S_2^2(\tau )}.$$

$$f_{\mathrm{SC}}=f_{\mathrm{PL}}+(2\pi {)}^{-1}\mathbf{k}\cdot {\mathbf{V}}_{\mathbf{rel}},$$

$$f_{\mathrm{e}}\approx \frac{\pi V_{\mathrm{SW}}^2}{V_{\mathrm{A}}(r_{\mathrm{SC}}-r_{\mathrm{source}})}.$$

Full text

\turnoffediting

1 Introduction

Voyagers 1 and 2 (V1 and V2) crossed the heliospheric termination shock (TS) in December 2004 and in August 2007, respectively (Stone et al., 2005, 2008). It is widely accepted that on 25 August 2012, V1 crossed the heliopause (HP) and is now moving through the local interstellar medium (LISM) (Burlaga et al., 2013b; Burlaga & Ness, 2014a; Gurnett et al., 2013; Stone et al., 2013; Webber & McDonald, 2013).

The magnetic field behavior beyond the HP is determined by the interstellar magnetic field (ISMF) draping over the HP as a tangential discontinuity that separates the solar wind (SW) from the LISM. The simulations of Pogorelov et al. (2009a, 2017b); Borovikov & Pogorelov (2014) predicted the magnetic field behavior very similar to observations. The recent simulation of Kim et al. (2017) shows that the ISMF “undraping” is consistent with the time-dependent processes occurring in the SW. It also reproduces major shocks propagating through the LISM that cause the plasma wave events observed by PWS between Nov 2012 to Jul 2016 (Gurnett et al., 2013, 2015; Pogorelov et al., 2017b). Interestingly, it also predicted a shock passing through V1 in Aug 2017, when plasma waves were observed by PWS. V2 remains in the inner heliosheath (IHS) measuring a velocity profile very different from that of V1 at the same distance from the Sun and gradually approaches the HP. Conventionally, the IHS is the SW region between the TS and the HP. The heliopause has “structure” clearly observed by V1. It was crossed in about 1 month, which gives width of about 0.3 AU. There are strong indications that this “structure” is due to the HP instability (Borovikov & Pogorelov, 2014) or magnetic reconnection (Schwadron & McComas, 2013; Pogorelov et al., 2017b), or both.

Theoretical and numerical studies of the SW–LISM interaction have a long history reviewed, e.g., in Zank (1999, 2015); Izmodenov et al. (2009); Izmodenov & Alexashov (2015); Opher (2016), and Pogorelov et al. (2017a). The TS is formed due to deceleration of the supersonic wind when it interacts with the HP and LISM counter-pressure. Modern models of the SW–LISM interaction take into account the effects of charge-exchange between ions and neutral atoms,coupling of the heliospheric magnetic field (HMF) and ISMF, and treat non-thermal (pickup) ions (PUIs) as a separate component (see Pogorelov et al., 2016, 2017b, and references therein).

As far as the model uncertainties are concerned, remarkably many observations of the SW/LISM bulk flow and average magnetic field have been reproduced by simulations. E.g., the deflection of the LISM neutral H atoms is on average in the $BV$-plane, which is defined by the LISM velocity and ISMF vectors, $\mathbf{V}_{\infty}$ and $\mathbf{B}_{\infty}$, in the unperturbed LISM (Izmodenov et al., 2005; Pogorelov et al., 2008, 2009b; Katushkina et al., 2015). Kinetic energetic neutral atom (ENA) flux simulations of Heerikhuisen et al. (2010, 2014); Heerikhuisen & Pogorelov (2011); Zirnstein et al. (2015a, b, 2016, 2017, 2018) reproduced the IBEX ENA ribbon using the $BV$-plane consistent with the hydrogen deflection plane. A number of simulations reproduce Voyager measurements. Pogorelov et al. (2017b) demonstrated that the distribution of density in the heliospheric boundary layer (a region of decreased plasma density on the LISM side of the HP) is in agreement with PWS data. These models also reproduce the H density at the TS derived from PUI measurements (Bzowski et al., 2009) and observed anisotropy in the 1–10 TeV galactic cosmic ray (GCR) flux (Schwadron et al., 2014; Zhang et al., 2014; Zhang & Pogorelov, 2016).

There are certain details in the observations and simulations that clearly demonstrate that turbulence bears an imprint of physical processes occurring in the IHS and LISM. These are related to the turbulent character of the SW both in front of the TS and in the IHS, the presence of the heliospheric current sheet (HCS) that separates the sectors of opposite HMF polarity, the variability of the boundary between the sector and unipolar HMF, and the possibility that instabilities and magnetic reconnection destroy the HCS, thus resulting in the HMF decrease. In fact, simulations imply that turbulence may be affecting the interaction pattern by facilitating magnetic reconnection and instabilities. Moreover, the LISM turbulence may be affected by shocks propagating through it.

This paper is an attempt to address these issues by performing a turbulence analysis of Voyager data. In particular, this study provides a spectral characterization of the various regimes of magnetic field fluctuations from the energy-injection range, through the inertial cascade, down to scales where kinetic effects start to affect the dynamics ($\approx 10^{4}$ km). The analysis cannot be further extended to smaller scales due to resolution- and accuracy-related issues of Voyager data. However, it is known that the inertial-cascade regime of turbulence keeps track of physical processes taking place at smaller scales, which makes its analysis significantly intriguing.

HMF sectors in the IHS. A characteristic property of the SW flow is the existence of a HCS that separates magnetic field lines of opposite polarity, which originate at the solar surface. The HCS propagates with the solar wind kinematically, provided that it has no back reaction on the flow. Theoretical and numerical, kinetic and multi-fluid analyses of magnetic reconnection across the HCS, have been discussed by Drake et al. (2010, 2017) and Pogorelov et al. (2013, 2017b). Magnetic reconnection may reveal itself as a tearing mode (or plasmoid) instability and may take place especially close to the HP, where the sector width decreases to negligible values.

Resolution of the sectors of alternating magnetic field polarity in the IHS is impossible. It is known that Voyager 1 had been observing negative radial velocity component for two years before it crossed the HP (Decker et al., 2012; Pogorelov et al., 2009a, 2012). Because of the piling up effect, the sector width should be negligible near the HP. Moreover, the sector width decreases to zero at the HCS tips, which makes attempts to resolve the traditional HCS structure very challenging. Besides, current sheets can be created not only due to the tilt of the Sun’s magnetic axis, but also due to stream interactions abundant in the heliosphere. The transition to chaotic behavior of the HMF occurs when the sector width becomes less than the numerical resolution (Pogorelov et al., 2017b). In nature, this is possible both due to pre-existing turbulence and magnetic reconnection across the current sheets, which also creates turbulence.

The approach followed by Borovikov et al. (2011) to track the HCS surface based on the assumption of HMF being unipolar is neither practical, nor acceptable. Numerical simulations allow us to determine what happens to B if the HMF is assumed to be unipolar: it is clear that the calculated magnetic field strength in this case is substantially overestimated as compared with V1 observations in the IHS (see the discussions, e.g., in Pogorelov et al., 2015, 2017b). Thus, the possibility of HMF depressions in the IHS covered by a sectored HMF should not be disregarded. It can be identified by the increased turbulence level in relevant regions. Voyager spacecraft provide us with appropriate measurements to answer these questions.

Richardson et al. (2016) have investigated the effect of the magnetic axis tilt on the number of HCS crossings and compared the observed and expected numbers. It has been shown that the number of HCS crossings substantially decreased two years after V1 and V2 entered the IHS. However, V2 might have entered the unipolar region at that time. It was concluded that there are indications of magnetic field decrease possibly due to magnetic reconnection across the HCS. In addition, as shown by Drake et al. (2017), V2 data reveal that fluctuations in the density and magnetic field strength are anticorrelated in the sector regions, as is expected from their magnetic reconnection modeling, but not in the unipolar regions. A possible dissipation of the HMF in such regions may also be an explanation of a sharp reduction in the number of sectors, as seen from the V1 data.

Turbulence in the heliosheath and LISM. An extensive data analysis related to the SW turbulence behavior in the SW ahead of the TS and in the IHS along the Voyager trajectories was performed by Burlaga (1994); Burlaga et al. (2003a, b, c, 2006a, 2007); Burlaga & Ness (2009); Burlaga et al. (2009, 2015, 2017, 2018). There is no single physical mechanism responsible for all observed turbulence manifestations. Large-amplitude fluctuations in the magnetic field strength $B$ are observed at small scales with very complex profiles. These fluctuations were described as “turbulence” (Burlaga et al., 2006a; Fisk & Gloeckler, 2008), although their nature and origin are not yet understood. The turbulence includes “kinetic scale” features (with sizes of the order of 10–100 gyroradii) and microscale features ($>100$ proton gyroradii). Usually, the observed turbulence consists of both coherent and random structures as seen in time profiles of the magnetic field strength on scales from 48 s to several hours. As shown by Burlaga et al. (2009), the large-scale (1 day) fluctuations measured at V1 in the unipolar and sector regions differ in some aspects: they have a log-normal distribution in the post-TS region, but Gaussian in the unipolar region. Instabilities and magnetic reconnection enhance turbulence, i.e, magnetic field dissipation in the sector region should naturally affect magnetic field and plasma density fluctuations inside the IHS.

Magnetic field fluctuations have been observed also in the LISM. However, those fluctuations are smaller than in the IHS (Burlaga et al., 2018) and their nature is not yet understood. In this paper, we investigate the LISM turbulence in different regions separated by shocks causing plasma wave emission observed by V1. Our methodology allows us to investigate these fluctuations in more detail than earlier. In fact, as discussed in §2, analysis of power spectra in the SW turbulence necessary to address these issues is a challenging task because of the sparsity of the 48 s data. After the TS, about 70% of magnetic field 48-s data are missing. Thus, sophisticated spectral estimation techniques become mandatory for obtaining reliable and physically meaningful results.

In the present study, Section 2 describes the data sets used for the analysis. Section 3 contains results of the analysis of the inner heliosheath, and it is split in two parts: §3.1 for *Voyager 2 * and §3.2 for Voyager 1. LISM turbulence is discussed in Section 4, and final remarks follow in Section 5. In Appendix A, we report information of the methods used for spectral estimation. Finally, Appendix B, contains information on variance anisotropy.

2 Voyager data in the IHS and LISM and selected periods.

This study considers different intervals in the inner heliosheath and in the local interstellar medium. In particular, we used magnetic field data at the highest resolution publicly available, the 48-s averaged data measured in situ by the Voyager Interstellar Mission (https://voyager.jpl.nasa.gov/mission/interstellar-mission/). In the heliosheath and beyond, the Voyager LMF magnetometers (MAG experiment, see Behannon et al. (1977)) sample magnetic field at a rate of 2.08 samples per second. The rate of the telemetry is 0.0208 Hz, and 48 s averages are periodically published in the NASA Space Physics Data Facility (https://spdf.gsfc.nasa.gov/) and can also be accessed via the COHO web site (https://cohoweb.gsfc.nasa.gov/coho/). Data are currently available through day-of-year (DOY) 271 of 2017 for V1, and through 2015 DOY 356 for V2, in RTN reference frame. The Heliographic RTN coordinate system is centered at the spacecraft. The R axis points radially outward from the Sun, the T axis is parallel to the solar equatorial plane and points in the direction of the Sun’s rotation, while the N axis completes the orthonormal triad. In the most sensitive LMF range, the level of noise is 0.006 nT. However, 1-sigma systematic errors due to the data calibration process and other sources of noise (sensors, electronics, telemetry system, ground tracking stations) are estimated around $\pm$0.02 and $\pm$0.03 nT at V1 and V2, respectively. Besides, the variability of errors makes such uncertainties rise up to $\pm 0.1$ nT in specific periods or for specific field components (Berdichevsky, 2009). In addition to the noise, limited telemetry coverage (the Camberra antennas of the CDSCC can only view V1 and V2 twelve hours per day) leads to data gaps of 8-16 hours per day. This point constitutes the major challenge for a spectral analysis, such as that presented for the first time in this study. The numerical techniques we used are synthetically described in Appendix A, and previously used in Fraternale et al. (2016); Gallana et al. (2016); Fraternale (2017).

In the IHS we consider four periods for V1 and for V2, respectively, after 2009 (see Figure 1a). In particular, at V1 we selected the intervals: (A1) 2009, DOY 22 - 2009, DOY 180 (109.5$\pm 0.77$ AU); (B1) 2010, DOY 180 - 2011, DOY 180 (115.65$\pm$1.79 AU); (C1) 2011, DOY 180 - 2011, DOY 276 (117.91$\pm$0.47 AU); (D1) 2011, DOY 276 - 2011, DOY 365 (118.81$\pm$0.44). During these periods, Voyager 1 sampled almost unipolar magnetic fields with northern “toward” polarity with respect to the Sun. Periods C1 and D1 are separated by a sector boundary crossing which occurred in 2011, DOY 276 and lasted about one day, when the polarity became southern, “away”, until 2012 DOY 209. During this period, interaction with the local interstellar plasma likely occurred. A detailed description of sector boundaries in proximity of the heliopause from 2011.5 can be found in Burlaga & Ness (2014b).

As the Plasma Science (PLS) instrument is not operational at V1, the bulk wind speed sometimes can be recovered using a Compton-Getting analysis from Low Energy Charged Particle experiment (LECP) and Cosmic Ray Subsystem (CRS) data (Krimigis et al., 2011). The heliosheath plasma has been provided by Richardson & Burlaga (2013); Richardson & Decker (2014), and Richardson & Decker (2015). For our analysis, it is particularly important to highlight that V1 (traveling at $V_{\mathrm{SC1}}\approx 17$ km s*-1*, 34.5∘ North) measured low radial velocity components since its crossing the heliospheric termination shock (TS). In particular, the radial velocity $V_{\mathrm{R}}$ decreased almost linearly from about 100 km s-1 (2006) to 0 km s-1 in 2010.5, while the tangential speed $V_{\mathrm{R}}$ oscillated around -40 km s-1. Numerical simulations (Pogorelov et al., 2013) suggest that the absence of sector boundary crossings observed by V1 near the heliopause (HP), as well as the negative radial velocity observed during 2011 could be symptom of V1 being inside a magnetic barrier.

The flow at Voyager 2 ($V_{\mathrm{SC2}}\approx 15$ km s-1, 30∘ South) was quite different. The bulk speed remained almost constant at about 150 km s-1 throughout the IHS. In contrast to V1 observations, high variability was found at V2 in the fluxes of energetic particles, as shown by Decker et al. (2008). Such variability has been related to the possibility for V2 to be alternatively inside the unipolar region (UHS) or inside the sector region (SHS) (Opher et al., 2011; Hill et al., 2014). The sector region is defined as the region swept by the heliospheric current sheet (HCS). Likely, V2 has been very close to the boundary between these two regions, according to models based on kinematic propagation of the HCS’s maximal extension, measured by the Wilcox Solar Observatory (http://wso.stanford.edu/). Richardson et al. (2016) compared the results of two models with the actual number of sector boundary crossings observed at V2. It is believed that the spacecraft remained in the unipolar region from until 2013.83 when it entered the sector region, as discussed by Burlaga et al. (2017).

Based on literature analysis, we selected four intervals for V2, see Figure 1(b): (SHS1) 2009, DOY 62 - 2009, DOY 210 (89.0$\pm$0.65 AU); (UHS1) 2010 DOY 252 - 2011 DOY 210 (94.08$\pm$1.89 AU); (UHS2) 2012, DOY 1 -2013, DOY 300 (100.20$\pm$2.88 AU); (SHS2) 2013, DOY 300 - 2015, DOY 1 (104.95$\pm$1.87 AU). The thermal ions average plasma density ($n_{i}$) was about 0.001 cm-3 from 2009 to 2012, than increased to 0.002 cm-3 until 2015 (standard deviation is 0.0007 in the last period and 0.0003 in the earlier periods). The thermal plasma temperature was 63500 K in SHS1, 44700 K in UHS1, 56200 K in UHS2, and 51700 K in SHS2 (standard deviation about 22000 K).

Ultimately, we considered four consecutive intervals of Voyager 1 data in the local interstellar medium (LISM), see Figure 1c): (L1) 2012, DOY 340 - 2013, DOY 130 (123.3$\pm$0.77 AU); (L2) 2013 DOY 133 - 2014 DOY 236 (126.4 $\pm$2.29 AU); (L3) 2014 DOY 273 - 2015 DOY 135 (130.2$\pm$ 1.11 AU); (L4) 2015 DOY 145 - 2016 DOY 246 (133.7$\pm$ 2.28 AU). This partition was also identified by Burlaga & Ness (2016): in that study, L1 and L3 were referred to as “disturbed” intervals while L2 and L4 as “quiet”. The periods L2 and L4 have been also consedered by Burlaga et al. (2015) and Burlaga et al. (2018), respectively. These regions are bounded by weak perpendicular shocks (or pressure waves) propagating through the LISM, as shown by Burlaga et al. (2013a); Gurnett et al. (2013, 2015); Burlaga & Ness (2016), and Kim et al. (2017). The interstellar plasma is colder than the IHS ($T\approx 10^{4}$ K). Electron plasma oscillations detected by V1’s Plasma Wave Subsystem (PWS) yielded the density estimate of 0.08 cm-3 (Gurnett et al., 2013). These oscillations are driven by electron beams produced upstream of the shocks. Until 2016, five events have been detected (Oct-Nov 2012, Apr-May 2013, Feb-Nov 2014, Sep-Nov 2015).

I this study, we removed outliers from the 48 s data sets. This was done by computing a backward and forward moving variance over a 48-hours window (3600 data points). For each magnetic field component, a data point was removed if it was larger than 6 times the minimum between the backward and the forward moving standard deviation ($|B_{j}|>6\min\{\sigma_{b_{j}},\sigma_{f_{j}}\},j=1,\dots,n$). Moreover, three calibration events (spacecraft rolls with 30 minutes periodicity) were removed from the intervals A1 and SHS1. Magnetic field data have been rotated to mean field coordinates: a B-parallel component $B_{\parallel}$, and two perpendicular components $B_{\perp 1}$, $B_{\perp 2}$ with respect to the average field $\mathbf{B_{0}}$, which better suits turbulence analysis . Since the SW flow and the magnetic field directions in the IHS are nearly orthogonal, $B_{\parallel}\approx B_{\mathrm{T}}$, $B_{\perp 1}\approx B_{\mathrm{N}}$ and $B_{\perp 2}\approx B_{\mathrm{R}}$. Moreover, for each component, we removed linear trends.

In the following, bold-shaped letters indicate vector fields and non-bold letters are used for the magnitude of vector fields.

3 Magnetic field fluctuations in the inner heliosheath

This section shows the results of magnetic field fluctuation analysis in the IHS. We provide a spectral analysis for a frequency range wider than five decades, $f\in[10^{-8},10^{-2}]$ Hz. Before discussing results, few definitions and symbols are introduced.

We computed the power spectral density ($\mathrm{PSD}$, or $P$), shown in Figures 2, 5 for the IHS, and in left panels of Figure 7 for the LISM. Due to the missing data issue, the PSD is estimated via three different numerical procedures (Appendix A). The comparative analysis of these techniques allows to recover the PSD with uncertainty of spectral indexes typically smaller 10%.

We investigated the spectral compressibility the ratio between the PSD of the field magnitude and the trace: $C(f)=\mathrm{P}[B]/E_{\mathrm{m}}$, where $E_{\mathrm{m}}(f)=\mathrm{tr}(\mathrm{P}[\mathbf{B}])=\mathrm{P}[B_{\parallel}]+\mathrm{P}[B_{\perp 1}]+\mathrm{P}[B_{\perp 2}]$. This measure is an index of the alignment of the fluctuation vector with the average field.

[TABLE]

In the above expressions, $\delta\mathbf{B}=\mathbf{B}-\mathbf{B_{0}}$ is the three-dimensional fluctuation about the background field, $\delta B=(\delta B_{\parallel}^{2}+\delta B_{\perp 1}^{2}+\delta B_{\perp 2}^{2})^{1/2}$ is its magnitude, $\mathbf{b}=\mathbf{B_{0}}/B_{0}$ is the direction cosines vector. Angle brackets indicate the ensemble average over all the data points of each interval, and the dot indicates the scalar product. The powers of 2 used in Eq. 1 allows us to interpret C1, C2 as the average percentage of fluctuating energy in the direction of $\mathbf{B_{0}}$. To reduce the contribution of noise, the computation is performed with hourly-averaged data. Compressibility values are reported in Tables 2, 4 (IHS) and 5 (LISM). It should be noted that for a fluctuating field with constant magnitude and isotropic angle distribution, C=0.33.

These tables also report information on other fluctuations properties such as the average turbulence intensity

[TABLE]

the maximum-variance fluctuation amplitude ($\delta B_{mv}$) and direction with respect to $B_{0}$ ($\theta_{mv}$), spectral breaks and spectral indexes, power-law exponents for the structure functions.

In fact, we performed a multi-scale analysis of the magnetic field increments via computation of the structure functions $S_{p}(f)$, a classical and powerful statistical tool to investigate the departure from self-similarity and the intermittent behavior of turbulence (Monin & Yaglom, 1971; Frisch, 1995; Politano & Pouquet, 1995; Politano et al., 1998).

[TABLE]

We used the absolute values in this definition for better convergence of statistics for the odd moments. The computation of $S_{p,j}$ from discrete data is nontrivial for Voyager data sets due to the amount and distribution of missing data. For this computation, we do not interpolate data and compute the statistics for the ensemble of available differences:

[TABLE]

where $\Delta t_{s}$ is the data resolution and $n$ the total number of points of the data set (we used both 48 s data and 1824 s averages). The counter $N(\tau_{j})$ of decreases with $\tau$ and also depends on the distribution of missing data. The amount of missing data is between 55% (A1) and 80% (D1), and the dominant periodicity of data gaps is 43000 s $\pm$ 2000 s ($f_{gap}=2.3\times 10^{-5}\pm 10^{-6}$ Hz). From the structure functions, one computes the scale-dependent kurtosis of magnetic field increments, which is an indicator of intermittency:

[TABLE]

Structure functions and kurtosis are shown in Figures 4 and 6 for the IHS, and 7 for the LISM.

The spectra computed from in situ single-spacecraft measurements are inevitably 1D-reduced spectra (Matthaeus & Goldstein, 1982) with frequencies measured in the spacecraft reference frame ($f_{\mathrm{SC}}$). The Doppler-shift relationship between the spacecraft and the plasma reference frame ($f_{\mathrm{PL}}$) reads

[TABLE]

where $\bm{k}$ is the vector wave number and $\mathbf{V_{rel}}=\mathbf{V_{SW}}-\mathbf{V_{SC}}$ the relative speed between the spacecraft and the plasma flow across it. Notice that for an Alfvénic nondispersive large-scale wave , the maximum value is reached for parallel fluctuations, $f_{\mathrm{PL}}=\kappa_{\parallel}V_{\mathrm{A}}/2\pi$. For dispersive waves instead (as in the kinetic regime), $f_{\mathrm{PL}}$ is typically a function of higher powers of the wavenumber. If $f_{\mathrm{PL}}\ll|\bm{\kappa}\bm{\cdot}\mathbf{V_{rel}}|/2\pi$, (Taylor’s hypothesis, Taylor (1938), see also the discussion in Howes et al. (2014)), frequencies measured at the spacecraft can be converted into wavenumbers in the direction of the relative wind flow. In the solar wind, this condition is satisfied on large scales as the flow is super-Alfvénic and the spacecraft is slow compared with the wind. This condition typically holds for the solar wind upstream of the TS, where $V_{\mathrm{SC}}/V_{\mathrm{SW}}\lesssim 0.05$. In the IHS, the situation differs quite at V1 and V2. At V2, in fact, $V_{\mathrm{SC}}/V_{\mathrm{SW}}\approx 0.1$ and $V_{\mathrm{A}}/V_{\mathrm{SW}}\approx 0.3$. The B-V angle is nearly equal to $\pi/2$, and the Taylor’s approximation might be used to obtain perpendicular wavenumbers, $\kappa_{\perp}\approx 2\pi f_{\mathrm{SC2}}/V_{\mathrm{SW}}$. We report the wavenumber value in all V2 figures. It should be reminded, however, that reduced spectra always contain contribution from all vector wavenumbers.

At V1, given that the spacecraft was in the slow-wind region, the Taylor’s approximation does not hold.

While performing the turbulence analysis of solar wind fluctuations, it is important to consider the causality condition. In fact, a fluctuation (or “eddy”) with typical velocity scale $\delta v$ and size $\ell$, experiences one “eddy turnover” in a period $t\sim\pi\ell/\delta v$. During this period, the eddy is convected by the wind by a distance equal to $d=V_{\mathrm{SW}}t$. Assuming the frozen-flow approximation, this fluctuation would be detected by the spacecraft at a frequency $f_{\mathrm{e}}\approx\kappa V_{\mathrm{SW}}/2\pi\approx V_{\mathrm{SW}}/\ell\approx\pi V_{\mathrm{SW}}/(\delta vt)\approx\pi V_{\mathrm{SW}}^{2}/(\delta v\ d)$. Reminding that $d$ is the distance traveled by the eddy from its origin to the spacecraft, one gets $f_{\mathrm{e}}\approx\pi V_{\mathrm{SW}}^{2}/[\delta v\ (r_{\mathrm{SC}}-r_{\mathrm{source}})]$. Eventually, considering $\delta v\approx V_{\mathrm{A}}$, the following expression is obtained

[TABLE]

This means that the frequencies less than $f_{e}$ in the spacecraft reference frame correspond to structures which did not experience yet one eddy turnover or, equivalently, waves that do not satisfy the causality condition, since they would have traveled from a further distance than their source point. In the solar wind upstream of the TS, the Sun can be clearly considered as the source point. In this case, fluctuations with $f_{\mathrm{SC}}\lesssim f_{e}$ would be older than the age of plasma. Fluctuations with $f_{\mathrm{SC}}\gtrsim f_{e}$ instead can be considered “active” fluctuations, meaning that they may be part of a turbulent energy cascade.

In this study, we computed $f_{e}$ at V2 only, and considered the termination shock as the source location.

Due to the differences highlighted above between the plasma flow at the two spacecraft and their spatial and temporal separation, magnetic field fluctuations in the IHS at V2 and V1 will be discussed in two separate subsections, §3.1 and §3.2, respectively.

3.1 IHS analysis of Voyager 2 data

Figure 2 shows magnetic field power spectra of the four selected periods of V2 data. Average plasma parameters are summarized in Table 1. Information about the average fluctuating energy, compressibility and strength of fluctuations is reported in Table 2, together with the frequency of spectral breaks and power-law exponents. Different regimes are identified.

Let us start the discussion with the high-frequency range ($10^{-3}\lesssim f<10^{-2}$ Hz). In principle, 48-s data could allow us to investigate the beginning of the kinetic regime, as the ion cyclotron frequency in the IHS is of the order of mHz (see Smith et al. (2006a); Alexandrova et al. (2008, 2012), and Schekochihin et al. (2009) for a detailed review). Unfortunately, however, noise affects the data, as explained in §2. Thus, all PSD figures contain a gray band at a power level of $P_{\mathrm{noise}}=0.029\ \mathrm{nT^{2}s}$, corresponding to a white noise of 0.03 nT amplitude (the actual distribution of the noise is unknown, this is a conservative threshold). The band is set at $P_{\mathrm{noise}}=0.086\ \mathrm{nT^{2}s}$ for $E_{\mathrm{m}}(f)$. The “noisy region” includes frequencies $f\gtrsim 5\times 10^{-4}$ Hz. Here, a spectral flattening in the last frequency decade the spectra show definite trends and retain some information on the anisotropy. Moreover, these trends are not strictly identical among data sets (see for instance the flattening and consecutive steepening at V1 during A1 in Figure 5(a)). .

Moreover, the cyclotron frequency of low-energy PUIs (1 keV) falls within this range, as shown in Table 1. It would be interesting to investigate the effect of PUIs in mediating turbulence and driving kinetic waves, which can affect the high-frequency part of the inertial regime, as shown in Smith et al. (2006b); Cannon et al. (2014b, a) and Aggarwal et al. (2016) for the solar wind from 1 to 6 AU. In fact, in the IHS the pickup-ion effect is expected to be considerable, since the density of 1-5 keV pickup protons is about 20% of the thermal protons density Zank et al. (2010).

Going back to larger scales, most spectra show the presence of a low-frequency regime ($f\lesssim 10^{-5}$ Hz) which we interpret as the Energy-Injection range (EI), a reservoir of energy for turbulence. Here, the magnetic energy decays as $1/f$. In particular, the spectral index $\alpha_{1}$ falls between 0.7 and 1.3 for all the components. We computed the spectral index by linear regression in the log-log space. The uncertainty due to the fit is always much smaller than that related to the different spectral estimation techniques. Errors are shown in Table 2, they are usually larger in the EI range.

Interestingly, a spectral break ($f_{b1}$) characterizes the end of the EI range and the beginning of a steeper cascade which will be interpreted as the Inertial Cascade (IC) of turbulence. The spectral break between the EI and the IC regimes and the large-scale $\sim 1/f$ power law are known to exist in the solar wind upstream the termination shock (first observation is by Tu et al. (1984)).

Interpretations for the nature of this regime include: the superposition of uncorrelated samples of solar surface turbulence having log-normal distributions of correlation lengths (Matthaeus & Goldstein, 1986) which can determine 1/f energy decay (Montroll & Shleisinger, 1982); the reflection of primarily outward-traveling Alfvén waves in presence of large-scale inhomogeneities (Velli et al., 1989; Perez & Chandran, 2013; Tenerani & Velli, 2017).

In the IHS, a $1/f$ energy decay was observed at V1 during 2009 by Burlaga & Ness (2012) via multi-fractal analysis. Their observations were limited to the range $f\in[10^{-7},10^{-5}]$ Hz, which did not allow them to investigate the existence and location of spectral breaks in the IHS, shown here for the first time (the only published power spectra in the IHS are shown in Burlaga & Ness (2010, 2012, 2014b)).

Figure 5 and Table 2 show that at V2 the break frequency $f_{b1}$ observed in $E_{m}(f)$ (black curves) is about $5\times 10^{-5}$ Hz for SHS1, corresponding to a spatial scale $\ell_{b1}\approx 0.2$ AU along the wind direction. Considering the age of fluctuations born at the Sun, one obtains from Eq. (7) a cutoff frequency of $10^{-8}$ Hz. The actual break is instead more consistent with $f_{e}$ of fluctuations generated at the TS (or affected by it). Moreover, it is not physically reasonable to consider turbulent structures with size greater than the outer scale of the system. The IHS width observed by V1 is around 27 AU, corresponding to $\kappa\approx 1.5\times 10^{-12}\ \mathrm{m}^{-1}$ (note the wavenumber axis $\kappa\sim\kappa_{\perp}$ in all V2 plots). Moreover, the Sun’s rotation acts as a forcing with $f_{\mathrm{sun}}\approx 4\times 10^{-7}$ Hz. The sector spacing should be around 2 AU after TS, decreasing as the HP is approached, even though the canonical sector structure is no longer recognizable beyond 10 or 20 AU.

It seems that $f_{b1}$ increases in the unipolar periods UHS1 and UHS2 to values around $10^{-5}$ Hz ($\ell_{b1}\approx 0.1$ AU), while it decreases again to $f_{b1}\approx 7\times 10^{-7}$ Hz ($\ell_{b1}\approx 1.5$ AU) in the sectored interval SHS2, where the break is actually very weak. Note also that the break location differs for the $\delta B_{\parallel}$, $\delta B_{\perp 1}$ and $\delta B_{\perp 2}$ components, represented by the red, green and blue curves, respectively, in Figure 2. In fact, the break of the B-perpendicular fluctuations occurs at a higher frequency with respect to the B-parallel one, by a factor between 2 and 7 for all periods except the last one.

The high slope in the power spectra at $f\gtrsim f_{b1}$ suggests that a turbulent inertial cascade is ongoing. Moreover, in all intervals except UHS1 a second spectral knee is observed at $f_{b2}\approx 10^{-4}$ Hz ($\ell_{b2}\approx 0.01$ AU). It is particularly visible in the $B_{\parallel}$ component, while it is weaker and not always observed in $B_{\perp}$. The spectral slope between the two breaks, $\alpha_{2}$, is about -1.6 for SHS1 and SHS2, close to the Kolmogorov’ s -5/3 value (Kolmogorov, 1941) and compatible with the model of Goldreich & Sridhar (1995, 1997) for the k-perpendicular cascade of critically balanced turbulence. The latter model, however, is not quite adequate, since it ignores the compressibility, which plays a fundamental role in heliosheath and interstellar turbulence.

During unipolar periods the magnetic energy decays slightly faster and the index is around -1.75. Note that, in general, the perpendicular components contribute to the steepening of the spectra, and that the UHS2 spectrum changes slope in rather a continuous way across frequencies. Beyond $f_{b2}$, the slope increases to $\alpha_{3}\approx-2$. The cyclotron frequency of thermal protons, $f_{\mathrm{cp}}$, is around 2 mHz (see Table 1). Moreover, structures with size of the Larmor radius convected across the spacecraft may affect the power spectra at spacecraft-frame frequencies close to $f_{\mathrm{cp\ SC}}\approx 10^{-2}$ Hz. Structures of size comparable with the ion inertial radius, similarly, should be detected at $f_{\mathrm{ip\ SC}}\approx 2.5\times 10^{-2}$ Hz. In all likelihood, the second break $f_{b2}$ is still within the MHD inertial range. However, we suggest that gyroradii of 1 keV PUIs may affect the turbulence at $f_{\mathrm{cp\ 1keV\ SC}}\approx 10^{-3}$ Hz (this frequency shifts to lower values as $E(eV)^{-1/2}$). We see that in the range $f_{b2}\lesssim f\lesssim 10^{-3}$ Hz a reduction of compressibility and intermittency takes place.

To simplify comparisons and prevent misunderstandings, we emphasize that our $P[B_{\perp 1}]$ is sometimes referred to as the “perpendicular spectrum”, whereas $P[B_{\perp 2}]$ corresponds to the “quasi-parallel” spectrum (Matthaeus et al., 1990; Bieber et al., 1996). The ratio of these two spectra is equal to 1 in case of pure slab turbulence and 1.67 in the case of pure 2D turbulence. The slab/2D model is not descriptive of the IHS or LISM turbulence, as it ignores compressible fluctuations.

The spectral compressibility is shown in Figure 3 (top panel),

The existence of a turbulent inertial range is further investigated by the analysis of structure functions of temporal increments of the magnetic field, $S_{p}(\tau)$, for $p={1,2,3,4}$ (Eq. 3, 4). Results are shown in Figure 4 in terms of B-parallel and B-perpendicular structure functions (left column) and kurtosis (right). Similarly to the neutral-fluid turbulence, under the assumption of homogeneity and isotropy, the inertial range of MHD turbulence is defined as the range of scales where the third-order longitudinal structure function displays a linear dependence for the Elsässer field, $S_{3}\sim\tau$. More generally, a linear dependence occurs for the longitudinal flux of energy, a result known as the Yaglom’s four-thirds law (Monin & Yaglom, 1971), extended to the MHD case by Politano & Pouquet (1998). This quantity provides information on the rate of dissipation of the turbulent energy and related plasma heating (Sorriso-Valvo et al., 2007; Carbone et al., 2009; Hadid et al., 2018; Sorriso-Valvo et al., 2018).

In the Kolmogorov’s description of (non-intermittent) isotropic and homogeneous turbulence in fluids, $S_{p}(\tau)\sim\tau^{p/3}$ (Kolmogorov, 1941). The presence of intermittency makes the actual scaling exponents of both the velocity and the magnetic fields deviate from the linear trend. For models and observations, the reader is referred to She & Leveque (1994); Politano & Pouquet (1995, 1998); Bruno & Carbone (2013), and Muller & Biskamp (2000). Table 1 in Politano et al. (1998) shows a comparison of measured exponents with different theoretical predictions.

The structure functions of magnetic field in the IHS at V2 display a power-law behavior in the range of time scales which approximately corresponds to the frequency range between the spectral breaks identified in Figure 2. In the left panels of Figure 4, the red curves stands for $S_{p,\parallel}$ and the blue curves stands for $(S_{p,\perp 1}+S_{p,\perp 2})/2$ (Eq. 3). The structure functions are computed from Eq. 4 and show oscillations related to data gaps. The counter $N(\tau)$ is indeed a oscillating function of $\tau$, and on average it decreases linearly with $\tau$ (Eq. 4). When it is less than a certain threshold (specifically, when $N(\tau)<0.25\max[N(\tau^{\prime})],\ \tau^{\prime}\in[\tau-48h,\ \tau+48h]$), the color of curves is switched to gray. These points were not used to compute the scaling exponents $\zeta_{p}$.

As seen from the line curvature, each scaling exponent changes continuously across scales. However, the energy-injection range is easily identified, as well as the effect of the solar rotation.

Relative exponents of parallel and perpendicular fluctuations are similar, and appear to be closer to the values typical of plasma velocity rather than magnetic field (Politano & Pouquet, 1995). These values are consistent with the presence of inertial-range intermittency, which is confirmed by the familiar profiles of kurtosis shown in the panels (e)-(h) of Figure 4. In fact, the intermittency increases with frequency and starts to increase at the beginning of the IC range, approximately at $f_{b1}$. The EI range is characterized by Gaussian values ($K\approx 3$), or even sub-Gaussian. In the inertial-cascade range, $K(\tau)$ rises up to 10. It seems that a damping of such growth can occur at some point within the IC regime, see e. g. the period SHS2. This is also observed at V1 (Figure 6, left panels). The evolution of spectral compressibility seen in Figure 3 (top), suggests the existence of a relationship between intermittency and compressibility, as has been shown by Alexandrova et al. (2008) in the kinetic regime. The decrease of compressibility and intermittency in the high-frequency range deserves further investigation, at present we could either interpret it as (i) an effect of data noise or (ii) physical reasons, as observed in Sorriso-Valvo et al. (2017) at 1 AU. Again, the effect of the pickup ion populations should be considered.

It is worth noticing that the intermittency of magnetic fluctuations in the heliosheath was investigated earlier in the framework of a multi-fractal formalism (Meneveau & Sreenivasan, 1987; Frisch, 1995), see e.g. Burlaga et al. (2006b); Macek et al. (2012); Burlaga & Ness (2010, 2013); Macek & Wawrzaszek (2013); Macek et al. (2014). Most of the published analyses, however, are focused on the magnetic field magnitude and consider scales larger than one day. Notable exceptions are presented by Burlaga & Ness (2009) and Burlaga & Ness (2013), who consider the probability distribution functions of the increments of magnetic field magnitude. Burlaga & Ness (2009) also used 48 s data and provided a description of different magnetic structures observed in the IHS. We point out that the low values of multi-fractal index reported, e. g., by Macek et al. (2014) should not be interpreted as a non-intermittent inertial range of turbulence. In fact, their range of scales corresponds to the EI range in the present study.

3.2 IHS analysis of Voyager 1 data

Figure 5 shows power spectra for the V1 intervals A1, B1, C1 and D1. Average quantities are reported in Table 3, while the fluctuation statistics, spectral breaks and slopes, and structure-function exponents are shown in Table 4. Spectral compressibility is shown in Figure 3 (bottom panel), and high-order statistics are shown in Figure 6.

At V1, the intensity of magnetic fluctuations is in general smaller compared to V2. The contribution to the fluctuating energy is largely due to the $\delta B_{\parallel}$ components, which results in much larger values of compressibility ($C_{1}\approx 0.6$), especially during 2009-2011.5.

Since the Taylor’s hypothesis does not hold at V1, we do not convert spacecraft frequencies to wavenumbers and all figures show the spacecraft-frame frequency axis only. In the early periods A1 and B1, the spectra show a marked difference between B-parallel and B-perpendicular fluctuations, $\delta B_{\parallel}$ being more energetic than $\delta B_{\perp}$ by a factor as high as five in the central decades of the spectrum. At lower frequencies, the spectral index of the total energy is $\alpha_{1}\approx-1.2$ (black curves in Figure 5), while beyond the break a fast steepening occurs, with the slopes as high as $\alpha_{2}\approx-2.3$. Again, the shape of the spectral trace is mainly due to $\delta B_{\parallel}$, which displays a rather fast cascade in the range $10^{-5}\lesssim f\lesssim 10^{-4}$ Hz ($\alpha_{2,\parallel}\approx-2.5$). The two perpendicular components behave similarly to each other, and on average they experience a Kolmogorov-like spectral decay in the whole range of frequencies, with index between -1.35 and -1.8. In later periods C1 and D1, which are closer to the HP boundary, the spectral break becomes weaker even for $\delta B_{\parallel}$, so that the discrepancy between the components is significantly reduced (see the bottom panels in Figure 5 ). In fact, the level of compressibility decreases to 0.4 during D1 at all frequencies.

At frequencies higher than about $10^{-4}$ Hz, the spectra flatten, due to the lower limit of data accuracy. The curved shape at $10^{-3}<f<10^{-2}$ Hz during A1, however, suggests that physical phenomena such as wave-particle interaction or PUI-driven turbulence may still be relevant in the signal.

The profiles of parallel structure functions (red curves in Figure 6, left column) show a well-defined change in the power law on the time scales corresponding to the observed spectral break. The exponents $\zeta_{p}$ reported in Table 4 have been fitted in the range $\tau\in[10^{4}\ s,\tau_{b1}]$ for A1 and B1, and $\tau\in[10^{4},2\times 10^{5}]$ s for C1 and D1. Thus, we show the exponents of $S_{p}$ for the $\alpha\approx-2$ part of the power spectra. Relative exponents instead hold for the whole range of frequencies. Also, at Voyager 1, the fluctuations are intermittent (right panels of Figure 6). In addition to previous large-scale analyses (e.g. Burlaga et al., 2006a; Macek et al., 2014), we show that the kurtosis profiles of magnetic increments increase with increasing frequency for all intervals. It should be noted, however, that the growth starts within the range $\tau\in[10^{5},10^{6}]$ s. This was expected for the $B_{\perp}$ components, but not so for $B_{\parallel}$: is seems that the intermittency of the compressible component starts prior to the energy spectral break, at about $f\approx 10^{-6}$ Hz, except for the interval A1.

4 LISM magnetic turbulence

Figure 7 shows the power spectra of the interstellar magnetic field in four V1 intervals L1–L4. . As shown in the summary Table 5, the intensity of magnetic fluctuations with respect to the background field ($B_{0\ \mathrm{LISM}}\approx 0.5$ nT) is nearly one order of magnitude smaller than in the IHS. The fluctuating magnetic energy increases significantly in the later periods L3 and L4, and L1 is the most quiet interval. There is little variation of the rms of parallel fluctuations among the intervals, even though it is higher during the central periods L2 and L3. In fact, the central periods are the most compressible ($C\approx 0.5$). We observe a progressive increase of transverse fluctuations in the $\perp_{2}$ (radial) direction. In fact, during L4 magnetic fluctuations are primarily transverse, especially at the largest scales, as highlighted in Figure 8 (pink curve) and in Figure 12(h). This fact was first pointed out by Burlaga et al. (2018).

Power spectra shown in the left panels of Figure 7 contain five frequency decades, a range unexplored so far. The figure shows the noise level corresponding to a white noise with 0.04 nT amplitude (gray bands), i. e. $P_{noise}\approx 0.05\ \mathrm{nT}^{2}$s. Burlaga et al. (2018) indicate that noise may affect the data at $f\gtrsim 4\times 10^{-5}$ Hz, which is consistent with the spectral flattening we observe. .

However, one can observe a spectral flattening (or a small bump) occurring, for all periods, in the range $10^{-6}<f<10^{-5}$ Hz. Moreover, the ion cyclotron frequency is smaller than in the IHS, $f_{ci,\mathrm{LISM}}\approx 10^{-4}$ Hz. At lower frequencies, the energy decays as a power law with spectral index close to the Kolmogorov’s -5/3 value. Values reported in Table 5 have been computed in the range $5\times 10^{-8}<f<3\times 10^{-6}$ Hz. As usual, errors indicate the discrepancy between the three spectral estimation techniques and it is higher than in the IHS, as only the first decade is considered. It should be noted that in the interval L2 a spectral flattening occurs at the lower frequencies ($f<3\times 10^{-7}$ Hz). This may indicate that the turbulence is young and locally generated or, alternatively, affected by local structures as shocks. During L4, magnetic fluctuations seem to change nature, since transverse fluctuations (both along $\perp_{1}$ and $\perp_{2}$) become dominant and a clear power-law decay of energy is observed, with spectral index about -1.9 which is in part due to some rapid shears in the signal. The parallel cascade is much slower, $\alpha\approx-1.4$.

The LISM turbulence is expected to show the features of intermittency, which has not been considered yet in the literature. Our results are shown in Figure 7 (e-h), which show the scale-dependent kurtosis of magnetic increments as was done for the IHS data in the previous section. We see that there is no intermittency at time scales $\tau\gtrsim 10^{6}$ s, for all intervals, since the kurtosis is smaller than three. At smaller scales, intermittency is observed for the $\perp_{1}$ component only, in the first three periods. During L4 instead, a significant increase of the kurtosis is observed for both $\perp_{1}$ and $\perp_{2}$, in the range $10^{5}\lesssim\tau\lesssim 3\times 10^{6}$ s, where the energy cascade is fast. It should be noticed that in the noisy range ($\tau\lesssim 10^{4}$ s) the statistic returns to Gaussian. Moreover, both during L1 and L4, reduction of $K$ starts at larger scales, at about $\tau\approx 3\times 10^{5}$ s. This corresponds to the flattening observed in the power spectrum. Though it is possible that such reduction is an artifact due to data uncertainty, it does not occur systematically, so at the current state of our analysis we do not exclude physical reasons. Increments of parallel fluctuations never show strong intermittency (differently to what was observed in the IHS).

The structure function exponents reported in Table 5 have been computed in the range $\tau=[5\times 10^{5},5\times 10^{6}]$ s. The absence of intermittency at larger scales may be a result of the passage of shocks - which could have made the pristine interstellar fluctuations more Gaussian - or a signature of locally produced turbulence. The possibility that MHD waves could be transmitted from the IHS to the LISM has also been hypothesized (Zank et al., 2017).

5 Summary and final remarks

This work provides a broadband spectral and high-order statistical analyses of magnetic field fluctuations in the IHS and LISM which give new insights into the properties of turbulence and its evolution. Our aim was to investigate the existence of different regimes of turbulent fluctuations, characterize them and describe their evolution in time and space.

We considered 12 data sets at different times and latitudes, after 2009, between 88 AU and 135 AU. Analysis of high-resolution (48 s) in situ measurements from both Voyager 1 and 2 spacecraft with the proposed advanced spectral-estimation techniques made it possible to investigate the evolution of fluctuations across more than five frequency decades ($10^{-8}<f<10^{-2}$ Hz), a range of scales which has not been explored in the literature so far. We focused our attention on the energy cascade, the compressible nature, anisotropy and intermittency of magnetic fluctuations.

In the IHS at V2, we identified the Energy-Injection and Inertial-Cascade regimes. Instead, the signatures of the kinetic regime were expected to show up in the last decade of the frequency spectrum ($10^{-3}<f<10^{-2}$ Hz), but at present .

The EI range is featured by a 1/f power law decay of magnetic energy, non-intermittent statistics of magnetic increments, and low compressibility. Its frequency extent depends on the observational interval considered, and seems to be wider in the unipolar periods. The first relevant scale highlighted in this study corresponds the EI/IC spectral break, and could not be just related to the nominal spacing of sectors ($\approx$2 AU), nor to the age of fluctuations generated or affected by the termination shock. In the unipolar periods, for instance, the spectral break occurs at $f\approx 10^{-5}$ Hz, corresponding to a spatial scale of about 0.15 AU along the wind direction. The originally more .

The IC regime is characterized by (i) steepening of power spectra towards values of spectral index close to the Kolmogorov’s value of -1.67, (ii) rapid growth of the kurtosis of magnetic field increments, and (iii) a clearly-defined power law decay of the third-order structure function $S_{3}$, with typical exponents of MHD turbulence. It can be concluded that unipolar periods showed a faster energy decay than sectored periods. The second typical scale in the V2 analysis was observed at $f\approx 10^{-4}$ Hz ($\ell\approx 5\times 10^{-3}$ AU). Here, a spectral knee occurs mainly for $\delta B_{\parallel}$, and the maximum anisotropy and compressibility is observed.

Fluctuations in the broad unipolar regions traveled by V1 before 2011.5 show a marked dominance of $\delta B_{\parallel}$. Parallel energy experiences a spectral break separating the large-scale regime with slow 1/f decay from the fast regime with $\alpha\approx-2.5$. Interpretation of V1 spectra is challenging due to the lack of accurate plasma velocity data and slow wind conditions, which do not allow us to compute the wavenumbers. It is possible that the parallel-wavenumber components due to Alfvénic or fast magnetosonic fluctuations dominate the spacecraft-frame frequency, but at present we cannot verify this hypothesis. We showed that V1 fluctuations are intermittent with a power-law increase of kurtosis anticipating a spectral break. A similar trend is also observed for the spectral compressibility, which reaches its maximum at the break frequency.

Finally, we have analyzed four LISM intervals. Even in this case, we used 48 s data to extend the range of frequencies considered in past literature studies and improve the accuracy of spectral estimates. We note that the level of compressible fluctuations is not higher than 0.6, and confirm the recently observed change of nature of turbulence during 2015/2016. Such transition consists mainly of an increase in perpendicular energy (especially, in the $\perp_{2}$ component). Moreover, in all intervals we observed a spectral flattening resembling a small bump, for $10^{-6}\lesssim f\lesssim 10^{-5}$ Hz. This bump also corresponds to an increase in compressibility. Its nature should be further investigated. Intermittency is mainly observed in transverse fluctuations. Such rapid increase begins at a spacecraft frequency of $3\times 10^{-7}$ Hz. It is then plausible that the observed LISM turbulence is locally modified by the periodic passage of shocks.

We expect our results to provide additional constraints on numerical and theoretical models of the outer heliosphere and, hopefully, shed light onto transport properties of energetic particles in these regions of space.

The authors acknowledge support from the MISTI-Seeds MITOR project “Spectral analysis of the solar wind beyond the termination shock - interstellar medium/heliosphere interactions” 2015-2016. FF was also supported by the post-doctoral grant “FOIFLUT” 37/17/F/AR-B. NP was supported, in part, by NASA grants NNX14AJ53G, NNX16AG83G, and 80NSSC18K1649NS, and NSF PRAC award OAC-1811176. JDR was supported under NASA contract 959203 from JPL to MIT. Computational resources for spectral and statistical analysis were provided by HPC@POLITO (http://www.hpc.polito.it).

Appendix A Methods for spectral analysis

Computing power spectra over a broad range of frequencies in the outer heliosphere is challenging because of the sparsity of the 48 s data (70% of magnetic field are missing). Fraternale et al. (2016), Gallana et al. (2016), Iovieno et al. (2016) and Fraternale (2017) have demonstrated that a careful application of different, independent, techniques makes it possible to recover the spectrum with proper accuracy (e.g., with the accuracy of 10% or lower in the spectral index.

1. Correlation method with linear data interpolation (CI): the PSD is obtained as the Fourier transform of the two-point auto-covariance function computed from linearly-interpolated data. It recovers the spectrum well in the low-frequency range, i. e. at frequencies smaller than the typical frequency of large data gaps ($f_{gap}\approx 2\times 10^{-5}$ Hz). Due to the low-pass effect of the linear interpolation, a spectral leakage is observed at higher frequencies and spectral exponents are typically overestimated .

2. Compressed sensing spectral estimation (CS): a recent paradigm we adopted from the signal processing and telecommunication area (Candes et al., 2006a, b; Donoho, 2006). This method does not interpolate data. It allows to recover exactly sparse signals (signals with few nonzero frequencies) even if fewer data points than those required by the Shannon’s principle are available. Testing CS on turbulent, gapped, data sets, we found it recovers well the spectrum, especially the high-frequency range ($f\gtrsim f_{\mathrm{gap}}$). Depending on the interval considered, it may lack accuracy in the neighborhood of $f_{gap}$. This typically shows up as a lack of energy around $f\in[3\times 10^{-6},10^{-5}]$ Hz (green curves in Figure 9). Also, a small peak around $f_{\mathrm{gap}}$ is sometimes observed (see the case of SHS2, Figure 2d). CS has also been recently exploited for the analysis of the spectrum of magnetospheric intervals out of a Kelvin-Helmholtz-instability event observed by the Magnetospheric Multiscale (MMS) mission in the Earth magnetosphere (Sorriso-Valvo et al., 2019).

3. Optimization method (OP): it is a simple algorithm which aims at minimizing errors in the CI analysis. This method is based on a genetic algorithm (Charbonneau, 1995, our code is based on the open-source code PIKAIA) which returns the piecewise-linear model spectrum, $P_{\mathrm{OP}}(f)$ as a result. Starting from an initial model spectrum, a synthetic signal is obtained by inverse Fourier transform. From this synthetic set, some data are removed according to the gap distribution in the Voyager data set. The power spectrum of this gapped synthetic set is successively computed by using the CI technique ($P_{\mathrm{sy,~{}CI}}$). The result is then compared to the CI spectrum of the Voyager data set ($P_{\mathrm{voy,~{}CI}}$). The difference between the two spectra ($F=\sum_{i}|P_{\mathrm{sy,~{}CI}}(f_{i})-P_{\mathrm{voy,~{}CI}}(f_{i})|$) is minimized by the optimization algorithm which, at each generation, modifies the control points in the model spectrum which is an approximation of the true spectrum of the Voyager signal (Fraternale, 2017, Chapt. 4). This method helps to estimate the error of the linear interpolation. Heuristically, OP proved to work well for statistically-homogeneous data sets representing physical phenomena with a continuous spectrum distributed over a broad range of scales. It cannot represent peaks in the spectrum but, if any, they are well identified by the other methods above (Figure 9, blue curves).

Power spectra displayed in Figures 2, 5, and 7 are build by using the average result of CI and CS for $f<10^{-5}$ Hz, and for $f>10^{-5}$ Hz only. The results of the OP method is also shown in all PSD pictures by a continuous smooth curve. In general, a very good agreement is observed. Note also that the CI and CS spectra are smoothed by means of a moving average with a constant-width window in the logarithmic space (i.e. the averaging points increase linearly with the frequency). This smoothing has no effect on the level of power and on the spectral index.

Figure 9 shows the application of methods CI (red curves), CS (green curves), OP (blue curves) to three data sets, where data have been artificially removed according to the same gap distribution of the period B1. The first data set is a synthetic turbulence set (obtained by inverse transform of a given power spectrum - black line - using random phases of the Fourier coefficients) with a constant index $\alpha=-5/3$ (Figure 9a). In the second case $\alpha=-1$, which represents a harder test case (panel b). In the third case, we used gap-free Ulysses data in the period 1990.82-1991.31. Panels (d)-(f) show the (smoothed) ratio between the estimated and the true spectrum $P/P_{\mathrm{true}}$.

It should be noticed that - due to the gap distribution of Voyager magnetic field time series in the IHS and LISM - the frequency range where spectral estimation is more critical is around $10^{-5}<f<10^{-4}$ Hz. The lower bound corresponds to $f_{\mathrm{gap}}\approx 2\times 10^{-2}$ Hz, the higher is linked to the length of contiguous subsets.

Appendix B Variance anisotropy