* [19] M. de León, D. Martín de Diego, and A. Santamaría-Merino.  Tulczyjew triples and Lagrangian submanifolds in classical field theories.  In Applied Differential Geometry and Mechanic, pages 21–47. Univ. of Gent, Academia Press, Gent, 2003.
  * [20] M. de León, J. Marín-Solano, and J. C. Marrero.  The constraint algorithm in the jet formalism.  Diff. Geom. Appl., 6(3):275–300, 1996.  10.1016/0926-2245(96)82423-5.
  * [21] M. de León, E. Merino, J. A. Oubiña, P. R. Rodrigues, and M. Salgado.  Hamiltonian systems on $k$-cosymplectic manifolds.  J. Math. Phys., 39(2):876, 1998.  10.1063/1.532358.
  * [22] M. de León, E. Merino, J. A. Oubiña, and M. Salgado.  Stable almost cotangent structures.  Bolletino Unione Mat. Ital. B (7), 11(3):509–529, 1997.
  * [23] M. de León, E. Merino, and M. Salgado.  $k$-cosymplectic manifolds and Lagrangian field theories.  J. Math. Phys., 42(5):2092, 2001.  10.1063/1.1360997.
  * [24] M. de León, I. Méndez, and M. Salgado.  $p$-almost tangent structures.  Rend. Circ. Mat. Palermo, 37(2):282–294, 1988.  10.1007/BF02844526.
  * [25] M. de León, I. Méndez, and M. Salgado.  Regular $p$-almost cotangent structures.  J. Korean Math. Soc., 25(2):273–287, 1988.  https://jkms.kms.or.kr/journal/view.html?spage=273&volume=25&number=2.
  * [26] M. de León and P. R. Rodrigues.  Methods of Differential Geometry in Analytical Mechanics, volume 158 of Mathematics Studies.  North-Holland, Amsterdam, 1989.  10.1016/s0304-0208(08)x7115-4.
  * [27] M. de León, M. Salgado, and S. Vilariño.  Methods of Differential Geometry in Classical Field Theories.  World Scientific, 2015.  10.1142/9693.
  * [28] A. Echeverría-Enríquez, A. Ibort, M. C. Muñoz-Lecanda, and N. Román-Roy.  Invariant forms and automorphisms of locally homogeneous multisymplectic manifolds.  J. Geom. Mech., 4(4):397–419, 2012.  10.3934/jgm.2012.4.397.
  * [29] M. Forger and L. G. Gomes.  Multisymplectic and polysymplectic structures on fiber bundles.  Rev. Math. Phys., 25(9):1350018, 2013.  10.1142/S0129055X13500189.
  * [30] M. Forger and S. Z. Yepes.  Lagrangian distributions and connections in multisymplectic and polysymplectic geometry.  Differ. Geom. Appl., 31(6):775–807, 2013.  10.1016/j.difgeo.2013.09.004.
  * [31] P. L. García.  The Poincaré–Cartan invariant in the calculus of variations.  In Symposia Mathematiea, volume 14, pages 219–246, London, 1974\. (Convegno di Geometria Simplettica e Fisica Matematica, INDAM, Rome, 1973), Academic Press.
  * [32] E. García-Toraño and T. Mestdag.  Conditions for symmetry reduction of polysymplectic and polycosymplectic structures.  2303.16037, 2023.
  * [33] H. Geiges.  An Introduction to Contact Topology, volume 109 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, 2008.  10.1017/CBO9780511611438.
  * [34] I. Gelfand, V. Retakh, and M. Shubin.  Fedosov manifolds.  Adv. Math., 104(1):104–140, 1998.  10.1006/aima.1998.1727.
  * [35] G. Giachetta, L. Mangiarotti, and G. A. Sardanashvily.  New Lagrangian and Hamiltonian Methods in Field Theory.  World Scientific, River Edge, 1997.  10.1142/2199.
  * [36] H. Goldschmidt and S. Sternberg.  The Hamilton-Cartan formalism in the calculus of variations.  Ann. Inst. Fourier, 23(1):203–267, 1973.  10.5802/aif.451.
  * [37] X. Gràcia, R. Martín, and N. Román-Roy.  Constraint algorithm for $k$-presymplectic Hamiltonian systems: application to singular field theories.  Int. J. Geom. Methods Mod. Phys., 6(5):851–872, 2009.  10.1142/S0219887809003795.
  * [38] X. Gràcia, X. Rivas, and N. Román-Roy.  Constraint algorithm for singular field theories in the $k$-cosymplectic framework.  J. Geom. Mech., 12:1–23, 2020.  10.3934/jgm.2020002.
  * [39] C. Günther.  The polysymplectic Hamiltonian formalism in field theory and calculus of variations I: The local case.  J. Diff. Geom., 25(1):23–53, 1987.  10.4310/jdg/1214440723.
  * [40] I. V. Kanatchikov.  Canonical structure of classical field theory in the polymomentum phase space.  Rep. Math. Phys., 41(1):49–90, 1998.  10.1016/S0034-4877(98)80182-1.
  * [41] J. Kijowski.  A finite-dimensional canonical formalism in the classical field theory.  Comm. Math. Phys., 30(2):99–128, 1973.  10.1007/BF01645975.
  * [42] J. Kijowski and W. M. Tulczyjew.  A symplectic framework for field theories, volume 107 of Lecture Notes in Physics.  Springer-Verlag, Berlin Heidelberg, 1st edition, 1979.  10.1007/3-540-09538-1.
  * [43] J. M. Lee.  Introduction to Smooth Manifolds, volume 218 of Graduate Texts in Mathematics.  Springer New York Heidelberg Dordrecht London, 2nd edition, 2013.  10.1007/978-1-4419-9982-5.
  * [44] P. Libermann and C.-M. Marle.  Symplectic Geometry and Analytical Mechanics, volume 35 of Mathematics and Its Applications.  Springer Dordrecht, 1987.  10.1007/978-94-009-3807-6.
  * [45] C.-M. Marle.  The Inception of Symplectic Geometry: the Works of Lagrange and Poisson During the Years 1808-1810.  Lett. Math. Phys., 90:3–21, 2009.  10.1007/s11005-009-0347-y.
  * [46] J. Marsden.  Darboux’s Theorem Fails for Weak Symplectic Forms.  Proc. Am. Math. Soc., 32(2):590–592, 1972.  10.2307/2037864.
  * [47] G. Martin.  A Darboux theorem for multi-symplectic manifolds.  Lett. Math. Phys., 16:133–138, 1988.  10.1007/BF00402020.
  * [48] M. McLean and L. K. Norris.  Covariant field theory on frame bundles of fibered manifolds.  J. Math. Phys., 41(10):6808–6823, 2000.  10.1063/1.1288797.
  * [49] E. E. Merino.  Geometría $k$-simpléctica y $k$-cosimpléctica. Aplicaciones a las teorías clásicas de campos.  PhD thesis, Universidade de Santiago de Compostela, 1997.
  * [50] J. Milnor.  On the existence of a connection with curvature zero.  Comm. Math. Helv., 32:215–223, 1958.  10.1007/BF02564579.
  * [51] L. K. Norris.  Generalized symplectic geometry on the frame bundle of a manifold.  In Proc. Symp. Pure Math., volume 54.2, pages 435–465. Amer. Math. Soc., Providence RI, 1993.  10.1090/pspum/054.2.
  * [52] A. M. Rey, N. Román-Roy, M. Salgado, and S. Vilariño.  On the $k$-symplectic, $k$-cosymplectic and multisymplectic formalisms of classical field theories.  J. Geom. Mech., 3(1):113–137, 2011.  10.3934/jgm.2011.3.113.
  * [53] X. Rivas.  Geometrical aspects of contact mechanical systems and field theories.  PhD thesis, Universitat Politècnica de Catalunya (UPC), 2021.  2204.11537.
  * [54] N. Román-Roy.  Multisymplectic Lagrangian and Hamiltonian Formalisms of classical field theories.  Symmetry Integr. Geom.: Methods Appl. (SIGMA), 5(100), 2009\.  10.3842/SIGMA.2009.100.
  * [55] L. Ryvkin and T. Wurzbacher.  An invitation to multisymplectic geometry.  J. Geom. Phys., 142:9–36, 2019.  10.1016/j.geomphys.2019.03.006.
  * [56] G. A. Sardanashvily.  Generalized Hamiltonian formalism for field theory. Constraint systems.  World Scientific Publishing Company, River Edge, NJ, 1995.  10.1142/2550.
  * [57] W. Ślebodziński.  Exterior forms and their applications.  PWN-Polish Scientific Publishers, Warsaw, 1970.  Revised translation from the French by Abraham Goetz. Monografie Matematyczne, Tom 52.
  * [58] S. Tchuiaga, F. Houenou, and P. Bikorimana.  On cosymplectic dynamics I.  Complex Manifolds, 9(1):114–137, 2022.  10.1515/coma-2021-0132.
  * [59] A. J. Tromba.  Almost-Riemannian structures on Banach manifolds: the Morse lemma and the Darboux theorem.  Can. J. Math., 28(3):640–652, 1976.  10.4153/CJM-1976-064-x.
  * [60] I. Vaisman.  Locally conformal symplectic manifolds.  Int. J. Math. Math. Sci., 8:521–536, 1985.  10.1155/S0161171285000564.
  * [61] I. Vaisman.  Connections under symplectic reduction.  Bull. Math. Soc. Sci. Math. Roum. N.S., 43(91):377–387, 2000\.  https://www.jstor.org/stable/43678676.


2022 290 20 December

Azimuthal anisotropy of jet particles in p–Pb and Pb–Pb collisions

ALICE Collaboration ALICE Collaboration

The azimuthal anisotropy of particles associated with jets (jet particles) at
midrapidity is measured for the first time in p–Pb and Pb–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}$ $=$ 5.02 TeV down to transverse momentum ($p_{\rm
T}$) of 0.5 GeV/$c$ and 2 GeV/$c$, respectively, with ALICE. The second-order
Fourier coefficient of the jet-particle azimuthal distribution ($v_{2}$) in
high-multiplicity p–Pb collisions is positive, with a significance reaching
6.8$\sigma$ at low $p_{\rm T}$. Comparisons with the inclusive charged-
particle $v_{2}$ and with AMPT calculations are discussed. The model describes
qualitatively the main features of the jet-particle $v_{2}$ in high-
multiplicity p–Pb collisions and indicates that the positive jet-particle
$v_{2}$ is generated by parton interactions.

The study of ultrarelativistic heavy-ion collisions aims to investigate the
properties of strongly-interacting matter characterised by high energy density
and temperature, known as the quark–gluon plasma (QGP) [1, 2]. In non-central
collisions, the initial spatial anisotropy of the overlap region of the
colliding nuclei is converted into an anisotropy in momentum space via
multiple interactions among the medium constituents. The final-state
anisotropies are quantified by the coefficients ($v_{n}$) of a Fourier
decomposition of the azimuthal ($\varphi$) distribution of produced particles
[3, 4]:

${{\rm d}^{2}N\over{{\rm d}p_{\rm T}{\rm d}\varphi}}={1\over 2\pi}{{\rm
d}N\over{\rm d}p_{\rm T}}\bigg{(}1+2\sum_{{n}=1}^{\infty}v_{n}(p_{\rm T}){\rm
cos}[{n}(\varphi-\Psi_{n})]\bigg{)},$ (1)

where $p_{\rm T}$ is the transverse momentum and $\Psi_{n}$ is the azimuthal
angle of the symmetry plane for the $n^{\rm th}$ harmonic. The azimuthal
anisotropy of produced particles is amongst the most important observables
used to probe the properties of the QGP. The second Fourier coefficient
$v_{2}$, referred to as elliptic flow [5, 4], provides information on the
collective expansion of the produced medium at low $p_{\rm T}$ [6] and the
path-length dependence of medium-induced parton energy loss at high $p_{\rm
T}$ [7, 8].

Collisions of small systems such as proton–nucleus are studied in detail to
characterise the cold nuclear matter (CNM) effects which are present in heavy-
ion collisions, providing critical information for understanding QGP
properties. These CNM effects influence parton distribution functions [9],
induce Cronin-like effects [10], and cause energy loss [11]. As in heavy-ion
collisions, a significant $v_{2}$ was observed in p–Pb collisions for soft as
well as hard probes, such as open heavy-flavour hadrons [12, 13], quarkonia
[14, 15], and high-$p_{\rm T}$ charged hadrons [16, 17]. A positive $v_{2}$
was also reported for charged particles [18, 19, 20, 21, 22] and muons from
charm-hadron decays [23] in pp collisions, while the $v_{2}$ of muons from
beauty-hadron decays was consistent with zero [23]. The observation of a non-
zero $v_{2}$ in small collision systems raised the question whether small-size
QGP droplets are formed in these conditions. However, the particle yields at
high $p_{\rm T}$ in p–Pb collisions are found to be unmodified compared to the
same measurements in scaled pp collisions within uncertainties, i.e. the
nuclear modification factor $Q_{\rm pPb}$ is compatible with unity [24, 25,
26, 27, 28, 17, 29]. Such an observation indicates that jet-quenching effects
are not significant in small collision systems. Therefore, alternative
scenarios have been proposed, see Ref. [30]. These scenarios include color
exchange in the final state [31], initial-state effects due to gluon
saturation [32, 30, 33], or anisotropic escape of partons from the surface of
the interaction zone [34].

This Letter reports the first measurements of the $p_{\rm T}$-differential
$v_{2}$ of jet particles (particles associated with jets) at midrapidity in
high-multiplicity p–Pb collisions and semicentral Pb–Pb collisions at a
centre-of-mass energy per nucleon pair $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with
ALICE. The jet-particle $v_{2}$ measurement in p–Pb collisions is of
particular interest since it extends down to lower $p_{\rm T}$ compared to
reconstructed jets and is fully separated from particles from soft processes.
These results are compared with inclusive charged-particle $v_{2}$
measurements in both p–Pb and Pb–Pb collisions, as well as with previous
results of the $v_{2}$ of reconstructed jets in Pb–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=2.76$ TeV [35]. Comparisons with AMPT predictions [36,
37] are also discussed.

The analysis is performed on the p–Pb and Pb–Pb data collected with the ALICE
detector in 2016 and 2015, respectively. In p–Pb collisions, the asymmetry of
the proton and Pb beam energies results in a rapidity shift of the
nucleon–nucleon centre-of-mass by 0.465 in the direction of the proton beam.
In the following, the pseudorapidity $\eta$ values correspond to the
laboratory frame. A detailed description of the detector and its performance
is given in Refs. [38, 39]. The analysis is based on tracks reconstructed
using the Time Projection Chamber (TPC) [40] located inside a large solenoidal
magnet with a 0.5 T field parallel to the LHC beam direction and covering
$|\eta|<0.9$. Information from the Inner Tracking System (ITS) [41] is used to
improve the spatial and momentum resolution of the reconstructed tracks. The
Silicon Pixel Detector (SPD) [42, 41], comprising the two innermost layers of
the ITS covering $|\eta|<2.0$ and $|\eta|<1.4$, is employed together with the
TPC to determine the position of the primary interaction vertex. The Forward
Multiplicity Detector (FMD) [43] consists of three sets of silicon strip
sensors, covering $-3.5<\eta<-1.7$ (FMD3) and $1.7<\eta<5$ (FMD1,2). The FMD
is used in p–Pb collisions for the event selection and to extract the $v_{2}$
coefficient via long-range three-particle correlations. The V0 detector,
formed by two scintillator arrays covering $-3.7<\eta<-1.7$ (V0C) and
$2.8<\eta<5.1$ (V0A), is used for triggering, event characterisation and
centrality determination [44]. Two sets of Zero Degree Calorimeters (ZDC)
[39], located at $\pm 112.5$ m from the interaction point along the beam line,
are also used for the event selection and characterisation.

The analysis of the p–Pb and Pb–Pb samples is based on events selected by a
minimum bias (MB) trigger. The MB trigger is provided by the coincidence of
signals in the two V0 scintillator arrays. Pile-up events are removed based on
an event selection which uses the information from the V0 and SPD to tag
events with multiple vertices. The beam-induced background is reduced offline
by exploiting the V0 and ZDC timing information. In addition, in p–Pb
collisions, the correlation between the multiplicity measured in the FMD and
V0 is used to further remove contamination from beam-induced background and
outliers in the FMD multiplicity distribution. Only events with a primary
vertex along the beam axis, $z_{\rm vtx}$, within $\pm 10$ cm from the nominal
interaction point are considered. About 526 million and 60 million p–Pb and
Pb–Pb events passed the event selection criteria, respectively. In Pb–Pb
collisions, the centrality classes are defined as percentiles of the Pb–Pb
hadronic cross section, and determined using the total energy deposited in the
V0 arrays [45]. The p–Pb data sample is divided into several multiplicity
classes based on the energy deposited in the V0A scintillators, located in the
Pb-going direction [45]. In the p–Pb analysis, the high-multiplicity class
0–10% and the low-multiplicity class 60–100% are studied. The measurements in
Pb–Pb collisions are presented in the 20–60% centrality interval.

Reconstructed charged-particle tracks are selected by applying the same
conditions on the number of space points, the quality of the track fit, and
the distance of closest approach to the primary vertex as in Ref. [46]. Tracks
are selected with $p_{\rm T}>0.5$ GeV/$c$ and $|\eta|<0.8$ in both p–Pb and
Pb–Pb collisions. Hits in the FMD are measured in the $\eta$ regions limited
to $-3.2<\eta<-1.8$ and $1.8<\eta<4.8$.

A detailed description of the procedure developed to calculate the jet-
particle $v_{2}$ is given in Ref. [47]. The jet-particle $v_{2}$ measurement
consists of four main parts: i) construction of the two-particle (charged
particles at midrapidity regarded as trigger and associated particles)
correlation function, ii) extraction of the near-side jet peak (the signal)
and background yields by fitting the two-particle correlation function, iii)
calculation of the $v_{2}$ of particle pairs represented by the trigger
particles using, for the first time, three-particle correlations in p–Pb
collisions and the scalar product method with the three-subevent technique in
Pb–Pb collisions, and iv) extraction of the $v_{2}$ of jet particles (signal),
in given trigger- and associated-particle $p_{\rm T}$ intervals, using a two-
component fit function that takes into account the relative contribution from
both jets and background to the particle-pair $v_{2}$.

The two-dimensional correlation function is constructed as a function of the
azimuthal angle difference ($\Delta\varphi$) and pseudorapidity difference
($\Delta\eta$) between trigger and associated particles (see Refs. [48, 14]
and references therein). Only pairs of particles with the same electric charge
are considered to suppress correlations originating from resonance decays. In
addition, the $p_{\rm T}$ of trigger particles ($p_{\rm T}^{\rm trig}$) is
taken to be larger than that of associated particles ($p_{\rm T}^{\rm assoc}$)
to avoid double counting of pairs. The correlation distribution is corrected
for the limited two-particle acceptance and detector inhomogeneities by using
the event-mixing technique [49], and normalized to the total number of trigger
particles.

Figure 1 (left) shows the two-dimensional correlation distribution in high-
multiplicity (0–10%) p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV
for $3<p_{\rm T}^{\rm trig}<5$ GeV/$c$ and $1<p_{\rm T}^{\rm assoc}<p_{\rm
T}^{\rm trig}$ GeV/$c$. A near-side jet structure is clearly observed at
$(\Delta\varphi\sim 0,\Delta\eta\sim 0)$ on top of the background. This
distribution is fitted with a double Gaussian and a sum of Fourier harmonics
up to the fifth order [50]. The former is used to extract the near-side jet-
peak yield ($S(\Delta\varphi,\Delta\eta)$), and the latter serves to obtain
the background yield ($B(\Delta\varphi,\Delta\eta)$).

Figure 1: Raw associated yield per trigger particle $Y$ (left) and $v_{2}$ of
trigger charged particles (right) as a function of $\Delta\eta$ and
$\Delta\varphi$ for pairs of charged particles measured in $|\eta|<0.8$ with
$3<p_{\rm T}^{\rm trig}<5$ GeV/$c$ and $1<p_{\rm T}^{\rm assoc}<p_{\rm T}^{\rm
trig}$ GeV/$c$ in high-multiplicity (0–10%) p–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV.

In the p–Pb analysis, the $v_{2}$ of particle pairs,
$v_{2}(\Delta\varphi,\Delta\eta$), is computed by using long-range three-
particle correlations. The particle pairs for a given
($\Delta\varphi,\Delta\eta$) cell are correlated with particles selected in
$1.8<\eta<4.8$ using FMD1,2 to construct the long-range correlation
distribution as a function of their azimuthal angle difference
($\Delta\varphi^{\prime}$) and pseudorapidity difference
($\Delta\eta^{\prime}$). Nonflow contributions, such as dijets, are suppressed
by subtracting the scaled ($\Delta\varphi^{\prime},\Delta\eta^{\prime}$)
correlation functions measured in low-multiplicity (60–100%) collisions
following the procedure described in Ref. [48]. In order to reduce the
statistical fluctuations of the subtracted
($\Delta\varphi^{\prime},\Delta\eta^{\prime}$) correlation functions, the
$\Delta\varphi^{\prime}$ projection is obtained from a first-order polynomial
fit along $\Delta\eta^{\prime}$ for each $\Delta\varphi^{\prime}$ interval
[51]. The correlation distribution in $-5.6<\Delta\eta^{\prime}<-1$ is further
projected onto $\Delta\varphi^{\prime}$ axis and fitted with a Fourier series
parameterised with the first three harmonics as

$\frac{{\rm d}N}{{\rm d}\Delta\varphi^{\prime}}\propto
1+2\sum_{n=1}^{3}V_{n\Delta}(\Delta\varphi,\Delta\eta){\rm
cos}(n\Delta\varphi^{\prime}),$ (2)

where $V_{2\Delta}(\Delta\varphi,\Delta\eta)$ is the relative second-order
harmonic. The procedure is repeated for each $p_{\rm T}$ interval of trigger
and associated charged particles. Assuming factorisation of single-particle
$v_{2}$ coefficients, the $V_{2\Delta}(\Delta\varphi,\Delta\eta)$ can be
expressed as the product of the $v_{2}$ of particle pairs
($v_{2}(\Delta\varphi,\Delta\eta)$) and the $v_{2}$ of particles in FMD1,2
($v_{2}^{\rm FMD1,2}$). The $v_{2}^{\rm FMD1,2}$ value, obtained with the
three-subevent technique [52] by constructing the long-range correlations
between the particles in the TPC, FMD1,2, and FMD3 detectors, amounts to
$0.028$ (with negligible uncertainties) in the 0–10% high-multiplicity class.
The $v_{2}(\Delta\varphi,\Delta\eta)$ distribution depicted in Fig. 1 (right
panel) shows a different structure in the region around the near-side jet peak
($\Delta\varphi\sim 0,\Delta\eta\sim 0$) compared to the background-dominated
region. This is a first indication of a different behaviour for
$v_{2}(\Delta\varphi,\Delta\eta)$ of jet and background particles.

In Pb–Pb collisions, the $v_{2}(\Delta\varphi,\Delta\eta)$ coefficient is
extracted from the scalar product method via the three-subevent technique [53,
54]. The method correlates particle pairs measured in the TPC with the second-
order event flow vector ${\bf Q}_{2}^{\rm{V0A}}$ estimated from the azimuthal
distribution of the energy deposited in the V0A. The resulting
$v_{2}(\Delta\varphi,\Delta\eta)$ is defined as

$v_{2}(\Delta\varphi,\Delta\eta)=\big{\langle}{\bf
u}_{2}(\Delta\varphi,\Delta\eta)\cdot{\bf
Q}_{2}^{\rm{V0A*}}\big{\rangle}_{\rm{TPC-TPC}}\bigg{/}\sqrt{\frac{\langle{\bf
Q}_{2}^{\rm{V0A*}}\cdot{\bf Q}_{2}^{\rm{V0C}}\rangle\langle{\bf
Q}_{2}^{\rm{V0A}}\cdot{\bf Q}_{2}^{\rm{SPD*}}\rangle}{\langle{\bf
Q}_{2}^{\rm{V0C}}\cdot{\bf Q}_{2}^{\rm{SPD*}}\rangle}},$ (3)

where ${\bf u}_{2}(\Delta\varphi,\Delta\eta)$ is the unit flow vector of each
particle pair measured in the TPC. The second-order harmonic event flow
vectors ${\bf Q}_{2}^{\rm{V0C}}$ and ${\bf Q}_{2}^{\rm{SPD}}$ measured in the
V0C and SPD, are introduced to take into account the resolution of the event
flow vector ${\bf Q}_{2}^{\rm{V0A}}$. The symbol * represents the complex
conjugate and the bracket $\langle...\rangle_{\rm{TPC-TPC}}$ denotes the
average over charged-particle pairs in a given $p_{\rm T}$ interval for
trigger and associated particles, and centrality range. The brackets in the
denominator indicate the average over all events in a centrality class
containing the particle pair. A recentering procedure is applied to correct
the event flow vectors for the non-uniform azimuthal acceptance effects of the
corresponding detector [55]. The pseudorapidity gaps between the TPC and V0A,
and the V0A, V0C, and SPD detectors suppress nonflow effects [56] and
eliminate autocorrelations [57].

In both p–Pb and Pb–Pb collisions, the $v_{2}(\Delta\varphi,\Delta\eta)$ can
be written [58, 59] as the weighted sum of the $v_{2}$ of jet particles
($v_{2}^{\rm jet\ part}$) and background ($v_{2}^{\rm B}(\Delta\varphi)$), as

$v_{2}(\Delta\varphi,\Delta\eta)=\frac{S(\Delta\varphi,\Delta\eta)}{S(\Delta\varphi,\Delta\eta)+B(\Delta\varphi,\Delta\eta)}v_{2}^{\rm
jet\
part}+\frac{B(\Delta\varphi,\Delta\eta)}{S(\Delta\varphi,\Delta\eta)+B(\Delta\varphi,\Delta\eta)}v_{2}^{\rm
B}(\Delta\varphi),$ (4)

where the jet-particle $S(\Delta\varphi,\Delta\eta)$ and background
$B(\Delta\varphi,\Delta\eta)$ yields are extracted from the two-particle
correlation functions constructed in the TPC. The $v_{2}^{\rm jet\ part}$
coefficient is obtained by parametrising $v_{2}^{\rm B}(\Delta\varphi)$ with a
Fourier series up to the fifth order and fitting Eq. 4 to the measured
$v_{2}(\Delta\varphi,\Delta\eta)$ distributions (Fig. 1, right), in a given
$p_{\rm T}$ interval for trigger and associated particles.

For comparison, the $v_{2}$ coefficient of inclusive charged particles is also
measured in p–Pb collisions using the three-subevent technique by constructing
long-range two-particle correlations, as done for the $v_{2}^{\rm FMD1,2}$
calculation, and in Pb–Pb collisions using the scalar product method.

The systematic uncertainties associated with the jet-particle $v_{2}$
measurement in high-multiplicity p–Pb collisions are estimated as follows. The
variation of the range of $z_{\rm vtx}$, which is decreased down to $\pm 8$
cm, and a less stringent condition on the correlation between the multiplicity
estimates obtained with the FMD and V0, give the systematic uncertainties
corresponding to the event selection. The uncertainty on the track
reconstruction is estimated by modifying the track selection criteria. The
bias due to the contribution of secondary particles produced in the FMD
acceptance on the jet-particle $v_{2}$ is investigated in A MultiPhase
Transport model (AMPT) simulations [36, 37]. In order to check for residual
nonflow effects after the subtraction of the scaled low-multiplicity event
class, the template fitting procedure [60] is tested. The difference between
the two procedures is considered as a systematic uncertainty. A potential bias
resulting from weak long-range correlations present in 60–100% low-
multiplicity events is studied by changing the interval from 60–100% to
70–100%. Finally, the $\Delta\varphi^{\prime}$ projection in Eq. 2 is obtained
from a constant fit instead of using a first-order polynomial fit along each
$\Delta\eta^{\prime}$ interval. The systematic uncertainties discussed above
are added in quadrature in each $p_{\rm T}$ interval of trigger and associated
particles to obtain a total systematic uncertainty on the jet-particle $v_{2}$
of 11.2–34.3%. The same systematic uncertainty sources affect the measurement
of the inclusive charged-particle $v_{2}$, resulting in a total systematic
uncertainty in the range 4.4–25.3%.

In Pb–Pb collisions, in addition to the systematic uncertainties arising from
the variation of the $z_{\rm vtx}$ range and track selection criteria listed
for p–Pb collisions, a systematic uncertainty related to the centrality
determination is estimated by using different centrality estimators for the
jet-particle $v_{2}$ measurement. The systematic effect related to the pile-up
event rejection is assessed via a dedicated analysis where pile-up events are
not removed, only to estimate their importance. The event flow vector is
computed by using V0C instead of V0A, giving a systematic effect on the
$v_{2}^{\rm jet\ part}$ calculation. All the systematic uncertainties are
added in quadrature to obtain the total systematic uncertainty ranging from
1.6–10.1%. The same sources of systematic uncertainties are associated to the
$v_{2}$ of inclusive charged particles and the total systematic uncertainty
varies up to 7.3%.

Figure 2 (top) presents the $v_{2}$ of jet particles as a function of $p_{\rm
T}$ at midrapidity ($|\eta|<0.8$) for different $p_{\rm T}^{\rm assoc}$
intervals, in high-multiplicity (0–10%) p–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV. The $p_{\rm T}$-differential inclusive
charged-particle $v_{2}$ coefficient is also displayed for comparison. A
positive jet-particle $v_{2}$ is observed with a significance of
2.6–6.8$\sigma$ for $p_{\rm T}\lesssim 5$ GeV/$c$, depending on the range of
$p_{\rm T}^{\rm assoc}$. The measured $v_{2}$ is independent of $p_{\rm T}$
within uncertainties and amounts to $\sim 0.04$. In contrast, the inclusive
charged-particle $v_{2}$ is larger in magnitude and shows a clear dependence
on $p_{\rm T}$. It increases up to $p_{\rm T}\sim 3$ GeV/$c$ where it reaches
a maximum value of $\sim 0.13$ and then decreases with increasing $p_{\rm T}$
to values similar to those of the jet-particle $v_{2}$.

Figure 2: Top: Jet-particle $v_{2}$ as a function of $p_{\rm T}$ of trigger
particles for several $p_{\rm T}^{\rm assoc}$ intervals in 0–10% p–Pb and
20–60% Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, compared with
the inclusive charged-particle $v_{2}$ in both systems. The values of the jet-
particle $v_{2}$ are horizontally shifted around the centre of the bin for
better visibility. The statistical uncertainties, shown as vertical bars, are
determined using the sub-sample technique. The systematic uncertainties are
represented as filled boxes. Horizontal bars indicate the bin width. The
previous results published by ALICE for reconstructed jet $v_{2}$ measured in
30–50% Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV are also shown
[35]. The Pb–Pb measurements are downscaled by a factor of 0.6. Middle:
Comparison of the jet-particle and inclusive charged-particle $v_{2}$ obtained
in p–Pb collisions with AMPT calculations [36, 37]. Bottom: Comparison of the
nuclear modification factor $Q_{\rm pPb}$ of inclusive charged particles
obtained in p–Pb collisions with AMPT calculations [36, 37].

A comparison with the jet-particle $v_{2}$ and the inclusive charged-particle
$v_{2}$ measured in semicentral (20–60%) Pb–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV is also reported in Fig. 2 (middle). The
Pb–Pb centrality class is chosen such that the eccentricity, which quantifies
the initial spatial anisotropies, is close to that in high-multiplicity p–Pb
collisions according to Glauber Monte Carlo simulations [61], although the
charged-particle multiplicities reached in the heavier system are larger. The
magnitude of the $v_{2}$ signal is also maximum in such collisions. The
inclusive charged-particle $v_{2}$ in Pb–Pb collisions is downscaled by an
empirical factor 0.6 in order to match the same measurement at $p_{\rm
T}\lesssim 3$ GeV/$c$ in p–Pb collisions. This factor, also applied to the
jet-particle $v_{2}$ and reconstructed-jet $v_{2}$ measured in Pb–Pb
collisions, may reflect the slightly different spatial anisotropies and the
larger multiplicities in Pb–Pb collisions [17]. In Pb–Pb collisions, the jet-
particle $v_{2}$ exhibits a mild decrease with increasing $p_{\rm T}$ from
$p_{\rm T}\sim 3$ GeV/$c$ and converges towards the $v_{2}$ of inclusive
charged particles for $p_{\rm T}\gtrsim 7$ GeV/$c$. The ALICE published
results for reconstructed jets measured in 30–50% Pb–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=2.76$ TeV [35] complement the present measurements by
extending them up to $p_{\rm T}=90$ GeV/$c$. They agree with recent
reconstructed-jet $v_{2}$ results obtained in Pb–Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV at $71<p_{\rm T}<298$ GeV/$c$ by ATLAS [62].
After this scaling, the jet-particle $v_{2}$ measured in p–Pb collisions has a
magnitude comparable to the $v_{2}$ of jet particles in the region $p_{\rm
T}\gtrsim 6$ GeV/$c$ and reconstructed jets measured in Pb–Pb collisions,
which are both attributed to jet-quenching effects [46]. In the intermediate
$p_{\rm T}$ region, $2<p_{\rm T}\lesssim 6$ GeV/$c$, the larger jet-particle
$v_{2}$ measured in Pb–Pb collisions is a consequence of interactions between
partons and the QGP medium.

The positive $v_{2}$ of jet particles measured in p–Pb collisions without any
indication of a modification of the jet production yields within experimental
uncertainties [17, 26, 25] may suggest that additional effects contribute to
the observed trends (Fig. 2, top). In order to shed more light on the origin
of the measured jet-particle $v_{2}$ in high-multiplicity p–Pb collisions,
Fig. 2 (middle) presents a comparison of the jet-particle and inclusive
charged-particle $v_{2}$ with the string-melting version of AMPT model [36,
37]. In this model, the $v_{2}$ is calculated following the same analysis
procedure as with data and the event characterisation is done by mimicking the
event class selection using the V0A detector at particle level, i.e. by
counting charged particles in $2.8<\eta<5.1$. The AMPT calculations generate a
positive jet-particle $v_{2}$ which shows a moderate dependence on the
associated-$p_{\rm T}$ interval. However, AMPT overestimates the measured jet-
particle $v_{2}$, predicting a $v_{2}$ whose shape and magnitude are
compatible with those of inclusive charged particles. Finally, the model
provides a fair agreement with the measured inclusive charged-particle
$v_{2}$. The enhanced jet-particle $v_{2}$ from AMPT calculations is
attributed to the missing string fragmentation process during the
hadronisation phase [63], in which only the coalescence mechanism is
incorporated [36]. As shown in Fig. 2 (bottom), the AMPT calculations provide
a nuclear modification factor $Q_{\rm pPb}$ of inclusive charged particles
[47] compatible with unity and consistent with the data at high $p_{\rm T}$
($p_{\rm T}\gtrsim 7$ GeV/$c$). From these comparisons, it can be concluded
that the azimuthal anisotropies in p–Pb collisions cannot be attributed to jet
quenching and are mainly driven by the non-equilibrium anisotropic parton
escape mechanism where partons have a higher probability to escape along the
shorter axis of the interaction zone [34]. In Pb–Pb collisions, the azimuthal
anisotropies are generated by the hydrodynamic expansion of the medium [64].

In summary, the jet-particle $v_{2}$ in high-multiplicity (0–10%) p–Pb
collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV is assessed for the first time
in the $p_{\rm T}$ range 0.5–8 GeV/$c$ by means of a novel multi-particle
correlation technique. The jet-particle $v_{2}$ is also measured down to low
$p_{\rm T}$ in semicentral Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02
TeV. These results are complementary to the previous jet $v_{2}$ results at
higher $p_{\rm T}$. A positive and $p_{\rm T}$-independent $v_{2}$ signal is
observed with a significance reaching 6.8$\sigma$ at low $p_{\rm T}$ in p–Pb
collisions. The $v_{2}$ magnitude is smaller than that measured in Pb–Pb
collisions at intermediate $p_{\rm T}$. The comparison with AMPT predictions
shows that parton interactions can generate a positive $v_{2}$ for jet
particles in high-multiplicity p–Pb collisions. These parton interactions
combined with the coalescence mechanism during the hadronisation phase, keep
the particle spectra unchanged with respect to binary-scaled pp collisions.
These new results bring crucial information about the origin of the observed
azimuthal anisotropies of jet particles in p–Pb collisions and set key
constraints on theoretical calculations.

## References

  * [1] HotQCD Collaboration, A. Bazavov et al., “Equation of state in (2+1)-flavor QCD”, Phys. Rev. D90 (2014) 094503, arXiv:1407.6387 [hep-lat].
  * [2] E. Shuryak, “Strongly coupled quark-gluon plasma in heavy ion collisions”, Rev. Mod. Phys. 89 (2017) 035001, arXiv:1412.8393 [hep-ph].
  * [3] S. Voloshin and Y. Zhang, “Flow study in relativistic nuclear collisions by Fourier expansion of Azimuthal particle distributions”, Z. Phys. C 70 (1996) 665–672, arXiv:hep-ph/9407282.
  * [4] A. M. Poskanzer and S. Voloshin, “Methods for analyzing anisotropic flow in relativistic nuclear collisions”, Phys. Rev. C 58 (1998) 1671–1678, arXiv:nucl-ex/9805001.
  * [5] J.-Y. Ollitrault, “Anisotropy as a signature of transverse collective flow”, Phys. Rev. D 46 (1992) 229–245.
  * [6] P. F. Kolb and U. W. Heinz, “Hydrodynamic description of ultrarelativistic heavy ion collisions”, SUNY-NTG-03-06 (2003) 634–714, arXiv:nucl-th/0305084.
  * [7] M. Gyulassy, I. Vitev, and X. Wang, “High $p_{\rm T}$ azimuthal asymmetry in noncentral A+A at RHIC”, Phys. Rev. Lett. 86 (2001) 2537–2540, arXiv:nucl-th/0012092.
  * [8] E. Shuryak, “Azimuthal asymmetry at large $p_{\rm T}$ seem to be too large for a pure ‘jet quenching”’, Phys. Rev. C 66 (2002) 027902, arXiv:nucl-th/0112042.
  * [9] K. J. Eskola, H. Paukkunen, and C. A. Salgado, “EPS09: A New Generation of NLO and LO Nuclear Parton Distribution Functions”, JHEP 04 (2009) 065, arXiv:0902.4154 [hep-ph].
  * [10] B. Z. Kopeliovich, J. Nemchik, A. Schafer, and A. V. Tarasov, “Cronin effect in hadron production off nuclei”, Phys. Rev. Lett. 88 (2002) 232303, arXiv:hep-ph/0201010 [hep-ph].
  * [11] Z.-B. Kang, I. Vitev, E. Wang, H. Xing, and C. Zhang, “Multiple scattering effects on heavy meson production in p+A collisions at backward rapidity”, Phys. Lett. B 740 (2015) 23–29, arXiv:1409.2494 [hep-ph].
  * [12] ALICE Collaboration, S. Acharya et al., “Azimuthal Anisotropy of Heavy-Flavor Decay Electrons in p–Pb Collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Rev. Lett. 122 (2019) 072301, arXiv:1805.04367 [nucl-ex].
  * [13] CMS Collaboration, A. M. Sirunyan et al., “Studies of charm and beauty hadron long-range correlations in pp and pPb collisions at LHC energies”, Phys. Lett. B 813 (2021) 136036, arXiv:2009.07065 [hep-ex].
  * [14] ALICE Collaboration, S. Acharya et al., “Search for collectivity with azimuthal J/$\psi$-hadron correlations in high multiplicity p–Pb collisions at $\sqrt{s_{\rm NN}}=$ 5.02 and 8.16 TeV”, Phys. Lett. B 780 (2018) 7–20, arXiv:1709.06807 [nucl-ex].
  * [15] CMS Collaboration, A. M. Sirunyan et al., “Observation of prompt J/$\psi$ meson elliptic flow in high-multiplicity pPb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 8.16 TeV”, Phys. Lett. B 791 (2019) 172–194, arXiv:1810.01473 [hep-ex].
  * [16] CMS Collaboration, A. M. Sirunyan et al., “Pseudorapidity and transverse momentum dependence of flow harmonics in pPb and PbPb collisions”, Phys. Rev. C 98 (2018) 044902, arXiv:1710.07864 [nucl-ex].
  * [17] ATLAS Collaboration, G. Aad et al., “Transverse momentum and process dependent azimuthal anisotropies in $\sqrt{s_{\mathrm{NN}}}=8.16$ TeV $p$+Pb collisions with the ATLAS detector”, Eur. Phys. J. C 80 (2020) 73, arXiv:1910.13978 [nucl-ex].
  * [18] ATLAS Collaboration, G. Aad et al., “Observation of Long-Range Elliptic Azimuthal Anisotropies in $\sqrt{s}=$13 and 2.76 TeV $pp$ Collisions with the ATLAS Detector”, Phys. Rev. Lett. 116 (2016) 172301, arXiv:1509.04776 [hep-ex].
  * [19] ATLAS Collaboration, M. Aaboud et al., “Measurements of long-range azimuthal anisotropies and associated Fourier coefficients for $pp$ collisions at $\sqrt{s}=5.02$ and $13$ TeV and $p$+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with the ATLAS detector”, Phys. Rev. C 96 (2017) 024908, arXiv:1609.06213 [nucl-ex].
  * [20] ATLAS Collaboration, M. Aaboud et al., “Measurement of multi-particle azimuthal correlations in $pp$, $p+$Pb and low-multiplicity Pb$+$Pb collisions with the ATLAS detector”, Eur. Phys. J. C 77 (2017) 428, arXiv:1705.04176 [hep-ex].
  * [21] CMS Collaboration, A. M. Sirunyan et al., “Observation of Correlated Azimuthal Anisotropy Fourier Harmonics in $pp$ and $p+Pb$ Collisions at the LHC”, Phys. Rev. Lett. 120 (2018) 092301, arXiv:1709.09189 [nucl-ex].
  * [22] ALICE Collaboration, S. Acharya et al., “Investigations of Anisotropic Flow Using Multiparticle Azimuthal Correlations in pp, p–Pb, Xe–Xe, and Pb–Pb Collisions at the LHC”, Phys. Rev. Lett. 123 (2019) 142301, arXiv:1903.01790 [nucl-ex].
  * [23] ATLAS Collaboration, G. Aad et al., “Measurement of azimuthal anisotropy of muons from charm and bottom hadrons in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”, Phys. Rev. Lett. 124 (2020) 082301, arXiv:1909.01650 [nucl-ex].
  * [24] ALICE Collaboration, J. Adam et al., “Centrality dependence of particle production in p–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Rev. C 91 (2015) 064905, arXiv:1412.6828 [nucl-ex].
  * [25] ALICE Collaboration, J. Adam et al., “Measurement of charged jet production cross sections and nuclear modification in p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV”, Phys. Lett. B 749 (2015) 68–81, arXiv:1503.00681 [nucl-ex].
  * [26] ALICE Collaboration, J. Adam et al., “Centrality dependence of charged jet production in p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV”, Eur. Phys. J. C 76 (2016) 271, arXiv:1603.03402 [nucl-ex].
  * [27] ATLAS Collaboration, G. Aad et al., “Transverse momentum, rapidity, and centrality dependence of inclusive charged-particle production in $\sqrt{s_{\rm NN}}=5.02$ TeV $p$ \+ Pb collisions measured by the ATLAS experiment”, Phys. Lett. B 763 (2016) 313–336, arXiv:1605.06436 [hep-ex].
  * [28] ALICE Collaboration, S. Acharya et al., “Measurement of electrons from heavy-flavour hadron decays as a function of multiplicity in p–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, JHEP 02 (2020) 077, arXiv:1910.14399 [nucl-ex].
  * [29] ALICE Collaboration, S. Acharya et al., “Measurement of prompt D0, D+, D∗+, and ${\mathrm{D}}_{\mathrm{S}}^{+}$ production in p–Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}=5.02$ TeV”, JHEP 12 (2019) 092, arXiv:1906.03425 [nucl-ex].
  * [30] K. Dusling, W. Li, and B. Schenke, “Novel collective phenomena in high-energy proton–proton and proton–nucleus collisions”, Int. J. Mod. Phys. E 25 (2016) 1630002, arXiv:1509.07939 [nucl-ex].
  * [31] K. Dusling and R. Venugopalan, “Evidence for BFKL and saturation dynamics from dihadron spectra at the LHC”, Phys. Rev. D 87 (2013) 051502, arXiv:1210.3890 [hep-ph].
  * [32] K. Dusling and R. Venugopalan, “Comparison of the color glass condensate to dihadron correlations in proton-proton and proton-nucleus collisions”, Phys. Rev. D 87 (2013) 094034, arXiv:1302.7018 [hep-ph].
  * [33] C. Zhang, C. Marquet, G.-Y. Qin, Y. Shi, L. Wang, S.-Y. Wei, and B.-W. Xiao, “Collectivity of heavy mesons in proton-nucleus collisions”, Phys. Rev. D 102 (2020) 034010, arXiv:2002.09878 [hep-ph].
  * [34] L. He, T. Edmonds, Z.-W. Lin, F. Liu, D. Molnar, and F. Wang, “Anisotropic parton escape is the dominant source of azimuthal anisotropy in transport models”, Phys. Lett. B 753 (2016) 506–510, arXiv:1502.05572 [nucl-th].
  * [35] ALICE Collaboration, J. Adam et al., “Azimuthal anisotropy of charged jet production in $\sqrt{s_{\rm NN}}=2.76$ TeV Pb–Pb collisions”, Phys. Lett. B 753 (2016) 511–525, arXiv:1509.07334 [nucl-ex].
  * [36] Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, “A Multi-Phase Transport Model for Relativistic Heavy Ion Collisions”, Phys. Rev. C 72 (2005) 064901, arXiv:nucl-th/0411110.
  * [37] H. Li, Z.-W. Lin, and F. Wang, “Charm quarks are more hydrodynamic than light quarks in final-state elliptic flow”, Phys. Rev. C 99 (2019) 044911, arXiv:1804.02681 [hep-ph].
  * [38] ALICE Collaboration, K. Aamodt et al., “The ALICE experiment at the CERN LHC”, JINST 3 (2008) S08002.
  * [39] ALICE Collaboration, B. Abelev et al., “Performance of the ALICE Experiment at the CERN LHC”, Int. J. Mod. Phys. A29 (2014) 1430044, arXiv:1402.4476 [nucl-ex].
  * [40] J. Alme et al., “The ALICE TPC, a large 3-dimensional tracking device with fast readout for ultra-high multiplicity events”, Nucl. Instrum. Meth. A 622 (2010) 316–367, arXiv:1001.1950 [physics.ins-det].
  * [41] ALICE Collaboration, K. Aamodt et al., “Alignment of the ALICE Inner Tracking System with cosmic-ray tracks”, JINST 5 (2010) P03003, arXiv:1001.0502 [physics.ins-det].
  * [42] ALICE Collaboration, ALICE Inner Tracking System (ITS): Technical Design Report.  Technical design report. ALICE. CERN, Geneva, 1999.  http://cds.cern.ch/record/391175.
  * [43] ALICE Collaboration, P. Cortese et al., ALICE technical design report on forward detectors: FMD, T0 and V0.  9, 2004.  https://cds.cern.ch/record/781854.
  * [44] ALICE Collaboration, E. Abbas et al., “Performance of the ALICE VZERO system”, JINST 8 (2013) P10016, arXiv:1306.3130 [nucl-ex].
  * [45] ALICE Collaboration, “Centrality determination in heavy ion collisions.” ALICE-PUBLIC-2018-011, Aug., 2018.  http://cds.cern.ch/record/2636623.
  * [46] ALICE Collaboration, S. Acharya et al., “Transverse momentum spectra and nuclear modification factors of charged particles in pp, p–Pb and Pb–Pb collisions at the LHC”, JHEP 11 (2018) 013, arXiv:1802.09145 [nucl-ex].
  * [47] ALICE Collaboration, “Supplemental material: Azimuthal anisotropy jet particles in p–Pb and Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, ALICE-PUBLIC-2022-020 (Dec, 2022) . https://cds.cern.ch/record/2845233.
  * [48] ALICE Collaboration, “Measurements of azimuthal anisotropies at forward and backward rapidity with muons in high-multiplicity p-Pb collisions at $\sqrt{s_{\rm NN}}=8.16$ TeV”, arXiv:2210.08980 [nucl-ex].
  * [49] ALICE Collaboration, B. Abelev et al., “Long-range angular correlations on the near and away side in p–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Lett. B 719 (2013) 29–41, arXiv:1212.2001 [nucl-ex].
  * [50] ALICE Collaboration, B. Abelev et al., “Multiplicity dependence of jet-like two-particle correlation structures in p–Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV”, Phys. Lett. B 741 (2015) 38–50, arXiv:1406.5463 [nucl-ex].
  * [51] ALICE Collaboration, J. Adam et al., “Forward-central two-particle correlations in p–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Lett. B 753 (2016) 126–139, arXiv:1506.08032 [nucl-ex].
  * [52] M. Luzum and J.-Y. Ollitrault, “Eliminating experimental bias in anisotropic-flow measurements of high-energy nuclear collisions”, Phys. Rev. C 87 (2013) 044907, arXiv:1209.2323 [nucl-ex].
  * [53] STAR Collaboration, C. Adler et al., “Elliptic flow from two and four particle correlations in Au+Au collisions at $\sqrt{s_{\rm NN}}$ = 130-GeV”, Phys. Rev. C 66 (2002) 034904, arXiv:nucl-ex/0206001.
  * [54] S. A. Voloshin, A. M. Poskanzer, and R. Snellings, “Collective phenomena in non-central nuclear collisions”, Landolt-Bornstein 23 (2010) 293–333, arXiv:0809.2949 [nucl-ex].
  * [55] I. Selyuzhenkov and S. Voloshin, “Effects of non-uniform acceptance in anisotropic flow measurement”, Phys. Rev. C 77 (2008) 034904, arXiv:0707.4672 [nucl-th].
  * [56] ALICE Collaboration, S. Acharya et al., “Energy dependence and fluctuations of anisotropic flow in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ and 2.76 TeV”, JHEP 07 (2018) 103, arXiv:1804.02944 [nucl-ex].
  * [57] ALICE Collaboration, S. Acharya et al., “Search for collectivity with azimuthal J/$\psi$-hadron correlations in high multiplicity p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 and 8.16 TeV”, Phys. Lett. B 780 (2018) 7–20, arXiv:1709.06807 [nucl-ex].
  * [58] ALICE Collaboration, S. Acharya et al., “Anisotropic flow of identified particles in Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, JHEP 09 (2018) 006, arXiv:1805.04390 [nucl-ex].
  * [59] ALICE Collaboration, E. Abbas et al., “J/$\psi$ Elliptic Flow in Pb-Pb Collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV”, Phys. Rev. Lett. 111 (2013) 162301, arXiv:1303.5880 [nucl-ex].
  * [60] ATLAS Collaboration, M. Aaboud et al., “Measurements of long-range azimuthal anisotropies and associated Fourier coefficients for $pp$ collisions at $\sqrt{s}=5.02$ and $13$ TeV and $p$+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with the ATLAS detector”, Phys. Rev. C 96 (2017) 024908, arXiv:1609.06213 [nucl-ex].
  * [61] C. Loizides, J. Kamin, and D. d’Enterria, “Improved Monte Carlo Glauber predictions at present and future nuclear colliders”, Phys. Rev. C 97 (2018) 054910, arXiv:1710.07098 [nucl-ex]. [Erratum: Phys.Rev.C 99, 019901 (2019)].
  * [62] ATLAS Collaboration, G. Aad et al., “Measurements of azimuthal anisotropies of jet production in Pb+Pb collisions at $\sqrt{s_{NN}}=$ 5.02 TeV with the ATLAS detector”, arXiv:2111.06606 [nucl-ex].
  * [63] W. Zhao, C. M. Ko, Y.-X. Liu, G.-Y. Qin, and H. Song, “Probing the Partonic Degrees of Freedom in High-Multiplicity $p-Pb$ collisions at $\sqrt{s_{NN}}$ = 5.02 TeV”, Phys. Rev. Lett. 125 (2020) 072301, arXiv:1911.00826 [nucl-th].
  * [64] A. Bzdak and G.-L. Ma, “Elliptic and triangular flow in $p$+Pb and peripheral Pb+Pb collisions from parton scatterings”, Phys. Rev. Lett. 113 (2014) 252301, arXiv:1406.2804 [hep-ph].


00 00 2012

# On the limitations of magneto-frictional relaxation

A. R. YEATES∗
Department of Mathematical Sciences ∗Corresponding author. Email:
<EMAIL_ADDRESS>Durham University  Durham  DH1 3LE  UK

###### Abstract

The magneto-frictional method is used in solar physics to compute both static
and quasi-static models of the Sun’s coronal magnetic field. Here, we examine
how accurately magneto-friction (without fluid pressure) is able to predict
the relaxed state in a one-dimensional test case containing two magnetic null
points. Firstly, we show that relaxation under the full ideal
magnetohydrodynamic equations in the presence of nulls leads necessarily to a
non-force-free state, which could not be reached exactly by magneto-friction.
Secondly, the magneto-frictional solutions are shown to lead to breakdown of
magnetic flux conservation, whether or not the friction coefficient is scaled
with magnetic field strength. When this coefficient is constant, flux is
initially conserved, but only until discontinuous current sheets form at the
null points. In the ensuing weak solution, we show that magnetic flux is
dissipated at these current sheets. The breakdown of flux conservation does
not occur for an alternative viscous relaxation scheme.

###### keywords:

Magnetohydrodynamics; magneto-friction; plasma relaxation

##  1 Introduction

Magneto-friction (hereafter abbreviated to MF) is a computational method for
obtaining force-free magnetic equilibria, whereby one retains the (ideal)
magnetohydrodynamic (MHD) induction equation

$\frac{\upartial\bm{B}}{\upartial
t}={\bm\nabla}\times\big{(}\bm{u}\times\bm{B}\big{)},$ (1)

but specifies the frictional velocity

$\bm{u}=\nu^{-1}\bm{J}\times\bm{B},$ (2)

rather than solving the full fluid equations. Here $\bm{B}$ is the magnetic
field and $\bm{J}=\mu_{0}^{-1}{\bm\nabla}\times\bm{B}$ is the current density.
This form of velocity changes equation (1) from hyperbolic to (degenerate)
parabolic type (Craig and Sneyd, 1986), and indeed the magnetic energy in a
volume $V$ satisfies

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{V}\frac{B^{2}}{2\mu_{0}}\,\mathrm{d}V=-\int_{V}\nu^{-1}|\bm{J}\times\bm{B}|^{2}\,\mathrm{d}V-\oint_{\upartial
V}\nu^{-1}B^{2}\bm{J}\times\bm{B}{\bm\cdot}\,\mathrm{d}\bm{S}.$ (3)

Thus energy is dissipated monotonically within the volume until a force-free
state $\bm{J}\times\bm{B}$ is reached.

The basic method dates back at least to Chodura and Schlüter (1981), and was
adopted (with $\nu=1$) by Craig and Sneyd (1986), who called it a “fictitious
fluid” and introduced a Lagrangian numerical method with the aim of better
preserving magnetic flux conservation. Yang et al. (1986) used a friction
coefficient of the form $\nu\sim B^{2}$, which has the effect of accelerating
the relaxation in regions of weak $B$. They observed a change in magnetic
topology during the relaxation, but attributed this to the effect of numerical
diffusion.

Subsequent authors have used MF in different ways to model the magnetic field
in the solar corona, taking it to be force-free to first approximation
(Wiegelmann and Sakurai, 2021). For example, Klimchuk and Sturrock (1992)
wrote $\bm{B}$ in terms of Euler potentials and used these to specify the
field line connectivity on the lower boundary (the solar photosphere).
Roumeliotis (1996) developed a “stress-and-relax” method for computing force-
free equilibria, where the magnetic field on the lower boundary is updated
step-by-step toward an observed vector magnetogram, with the coronal field
responding step-by-step by MF. This technique has been successfully used for
modelling solar active regions (Valori et al., 2005, 2010; Guo et al., 2016).

MF has also been used (with sufficiently small $\nu$) to model the quasi-
static response of the coronal magnetic field to slowly evolving boundary
conditions, rather than computing a single equilibrium. This idea was
introduced by van Ballegooijen et al. (2000), who modelled the evolution of
the large-scale coronal magnetic field in response to both large-scale solar
surface motions (primarily differential rotation) and small-scale
(supergranular) convection, which they parametrised by a surface diffusion
term. They showed that this simple model can capture the formation of sheared
filament channels, and even the sudden eruption of twisted magnetic structures
due to loss of equilibrium (e.g., Mackay and Yeates, 2012; Lowder and Yeates,
2017; Bhowmik and Yeates, 2021). The MF model has also been applied with
resolved small-scale boundary driving, using either an imposed convective
velocity field (Meyer and Mackay, 2016) or an imposed electric field
constrained by observed magnetograms (Mackay et al., 2011; Cheung and DeRosa,
2012; Pomoell et al., 2019; Hoeksema et al., 2020; Yardley et al., 2021). The
electric-field driven MF simulations have been successful at reproducing the
formation of observed non-potential magnetic structures within individual
active regions, and have also performed favourably in reconstructing the
injected magnetic energy and relative helicity when compared directly to a
full-MHD flux emergence simulation (Toriumi et al., 2020).

In this paper, our aim is to explore how well MF is able to reproduce the
equilibrium state that would be predicted by relaxation under the full ideal
MHD equations. Goldstraw et al. (2018) have recently published just such a
test for a setup that models a solar coronal loop whose magnetic field lines
are sheared by footpoint displacements at either end. They find that MF gives
an excellent approximation of the quasi-static sequence of MHD equilibria for
small plasma-beta, defined as $\beta=2\mu_{0}p/B^{2}$, where $p$ is the fluid
pressure. For large $\beta$, their MF result departs from the MHD solution
because it neglects the effects of pressure. Here, we build on that study by
considering a magnetic field containing null points ($B=0$). Near these
points, which are common in the coronal models mentioned above, the $\beta$
parameter is necessarily large and this raises the question of whether MF – as
given by equation (2) – will be successful at reproducing the ideal MHD
relaxation of such magnetic fields.

It is important to note that some implementations of the MF method do include
an additional fluid pressure term in equation (2), writing
$\bm{u}=\nu^{-1}(\bm{J}\times\bm{B}-{\bm\nabla}p)$. These include the original
calculations of Chodura and Schlüter (1981), those of Linardatos (1993), and
the Lagrangian schemes of Craig and Litvinenko (2005) or Candelaresi et al.
(2015). There are also implementations that add an MF term to a momentum
equation where the inertial terms are retained (Hesse and Birn, 1993;
Candelaresi et al., 2015). As such, our conclusions will apply only to the MF
model without additional pressure or inertial terms, which is nevertheless
widely used in solar physics.

Some particular criticisms of the MF method were raised by Low (2013). The
first important objection is that null points cannot move during magneto-
friction, since at such points equation (2) implies that $\bm{u}=\bm{0}$. This
suggests that the method will fail if the null points need to move during the
relaxation, as will be the case for our test in this paper. The second
objection is that discontinuous current sheets (jumps in $\bm{B}$ with
infinite $\bm{J}$) will form in finite time, so that the method breaks down.
We shall return to both of these important points below. Again, it should be
noted that both of these criticisms apply to the MF model without plasma
pressure, but need no longer apply if a pressure term is added to the
velocity.

The outline of the paper is as follows. Section 2 describes our test setup,
along with the “ground truth” ideal MHD solution. In order to be sure that our
results are not affected by numerical diffusion, we adopt a simple one-
dimensional setup, inspired by Bajer and Moffatt (2013). In section 3, we
present one-dimensional MF results for different forms of $\nu$, and we
conclude in section 4 with a discussion of the implications for future
modelling.

##  2 Test setup

Figure 1: The initial magnetic field for our test setup (dark line). The feint
line shows the state of minimum magnetic energy, discussed in section 2.2. The
positions of the null points in the former and current sheets in the latter
are indicated.

We will consider the relaxation of a magnetic field that initially has the
(non-equilibrium) form

$B_{x}(x,0)=B_{z}(x,0)=0,\hskip 28.45274ptB_{y}(x,0)=\tfrac{1}{2}+\sin(\pi
x),$ (4e,f)

on the domain $-1\leq x\leq 1$. For simplicity we will assume periodic
boundary conditions on all variables. This initial magnetic field is shown by
the solid curve in figure 1. The offset of $1/2$ is added so that the magnetic
null points, which are initially at $x=x_{0}=-1/6$ and
$x=x_{0}^{\prime}=-5/6$, will move during the relaxation, allowing a fairer
test of the MF method.

###  2.1 Ground truth: full MHD

Here we illustrate the solution of the full ideal MHD equations with our
initial condition (1), where we take the initial pressure and density to be
constant, with $p(x,0)=0$ and $\rho(x,0)=10^{-2}$. Thus we begin with a low-
beta plasma away from the locations where $B_{y}\approx 0$, consistent with
the solar coronal plasmas where magneto-friction is typically applied. We
assume an adiabatic ideal gas with polytropic index $\gamma=5/3$, so that the
ideal MHD equations for this initial condition reduce to

$\displaystyle\frac{\upartial{\rho}}{\upartial{t}}$
$\displaystyle=-\frac{\upartial{}}{\upartial{x}}\big{(}\rho u_{x}\big{)},$
(5a)
$\displaystyle\rho\left(\frac{\upartial{u_{x}}}{\upartial{t}}+u_{x}\frac{\upartial{u_{x}}}{\upartial{x}}\right)$
$\displaystyle=-\frac{\upartial}{\upartial{x}}\left(p+\frac{B_{y}^{2}}{2\mu_{0}}\right)+\mu\frac{\upartial{{}^{2}u_{x}}}{\upartial{x^{2}}},$
(5b)
$\displaystyle\frac{\upartial{p}}{\upartial{t}}+u_{x}\frac{\upartial{p}}{\upartial{x}}$
$\displaystyle=-\gamma
p\frac{\upartial{u_{x}}}{\upartial{x}}+(\gamma-1)\mu\left(\frac{\upartial{u_{x}}}{\upartial{x}}\right)^{\\!2},$
(5c) $\displaystyle\frac{\upartial{B_{y}}}{\upartial{t}}$
$\displaystyle=-\frac{\upartial{}}{\upartial{x}}\big{(}u_{x}B_{y}\big{)}.$
(5d)

A uniform viscosity $\mu$ is included, together with corresponding viscous
heating term in the pressure equation (5c) to ensure conservation of total
energy $E=E_{\rm kin}+E_{\rm int}+E_{\rm mag}$, where the kinetic, internal,
and magnetic energies are

$E_{\rm kin}=\int_{-1}^{1}\frac{\rho u_{x}^{2}}{2}\,\mathrm{d}x,\hskip
28.45274ptE_{\rm int}=\int_{-1}^{1}\frac{p}{\gamma-1}\,\mathrm{d}x,\hskip
28.45274ptE_{\rm
mag}=\int_{-1}^{1}\frac{B_{y}^{2}}{2\mu_{0}}\,\mathrm{d}x.\hskip 28.45274pt$
(6g–i)

We solve equations (5) numerically with the ATHENA code (Stone et al., 2008),
using the HLLD Riemann solver with third-order reconstruction on a uniform
grid of 1024 points in the $x$-direction. We set $\mu_{0}=1$ so that the
maximum Alfvén speed is order $\sqrt{\rho^{-1}}=10$ in our units. Results are
illustrated in figures 2(a,b) for the inviscid case $\mu=0$ as well as three
increasing values of viscosity $\mu$. Figures 2(c,d) show the full time
evolution for the single case $\mu=10^{-1}$.

Figure 2: Illustration of the full MHD relaxation. Panels (a) and (b) show
$B_{y}$ and $p$ at $t=10$ for several values of $\mu$, while panels (c) and
(d) show the time evolution up to $t=4$ of $B_{y}$ and $p$ for the case
$\mu=10^{-1}$. The thick gray line in (a) is the minimum energy solution in
section 2.2. (Colour online)

First, consider the evolution of $\bm{B}$. When the evolution begins, Alfvén
waves are immediately launched in opposite directions from the region of peak
$B_{y}(x,0)$ around $x=1/2$. These are seen as diagonal fronts in figure 2(c),
moving at speed approximately 10. These waves interact each time they cross
the domain, as the magnetic field relaxes to a lower energy state.

Figure 3: Time evolution of the three components of the total energy in the
full MHD relaxation, for different values of viscosity $\mu$. (For each
component, lighter coloured curves indicate larger $\mu$.) The thick gray
horizontal line (just visible) indicates the minimum magnetic energy under an
ideal relaxation with no further constraints (section 2.2). (Colour online)

Figure 3 shows that the magnetic energy reduces from its initial value of
$0.75$ to approximately $0.52$ during this initial relaxation phase, which
lasts until around $t=1$. After this time, the kinetic energy is very small,
even for $\mu=0$ although in that case some residual waves remain, as
evidenced by the oscillations at the Alfvén frequency in $E_{\rm mag}$. These
oscillations are damped by viscosity, and for $\mu=10^{-1}$ this happens
before $t=4$, as seen in figure 2(c).

It is clear from figure 3 that the magnetic energy lost during the relaxation
goes entirely into internal energy, as must be the case if a stationary
equilibrium is reached with $E_{\rm kin}=0$, owing to conservation of total
energy. This internal energy is manifested in the relaxed state by a non-zero
fluid pressure, as shown in figure 2(b). The pressure is strongest in the
regions where $B_{y}$ changes rapidly, consistent with this being a
magnetostatic equilibrium, which in our one-dimensional system amounts to a
total pressure balance

$\frac{\upartial{}}{\upartial{x}}\left(p+\frac{B_{y}^{2}}{2\mu_{0}}\right)=0.$
(7)

Thus in the final equilibrium the plasma-beta is quite variable: it is small
in the regions where $B_{y}$ was initially strong and so decreased, while it
is larger (0.2 to 0.4) in the region where the fluid has been compressed in
order to increase $B_{y}$. This conclusion is little changed by the presence
of viscosity (at least for the $\mu$ values considered here).

###  2.2 State of minimum magnetic energy

Since the MF method will be ignoring the fluid equations (5a-c), it is useful
to compare the full MHD solution to the state of minimum $E_{\rm mag}$
constrained only by the induction equation (5d). This would be a state of
uniform magnetic pressure, $B_{y}^{2}=\textrm{constant}$, and for our initial
condition (1), conservation of total unsigned flux shows that

$|B_{y}|=\frac{1}{2}\int_{-1}^{1}|B_{y}(x,0)|\,\mathrm{d}x=\frac{1}{6}+\frac{\sqrt{3}}{\pi}.$
(8)

Since the magnetic topology is preserved, the two null points must survive and
become discontinuous current sheets in the relaxed state, at which $B_{y}$
changes sign (cf. Bajer and Moffatt, 2013). From conservation of flux between
the nulls, and symmetry, it follows that these current sheets would be located
at

$x_{0}(T)=\frac{-1+\varDelta}{2},\hskip
28.45274ptx_{0}^{\prime}(T)=\frac{-1-\varDelta}{2},\hskip
28.45274pt\varDelta=\frac{6\sqrt{3}-2\pi}{6\sqrt{3}+\pi}.\hskip 14.22636pt$
(9g–i)

This minimum $E_{\rm mag}$ magnetic field is shown by the feint line in figure
1, repeated in figure 2(a) for comparison to the full MHD solution. It is
evident that the finite pressure in the full MHD solution modifies the shape
of the equilibrium and smooths out the current sheets to some extent.
Nevertheless, the magnetic energy in the MHD solution reaches close to that of
the minimum $E_{\rm mag}$ solution, which is

$\overline{E}_{\rm
mag}=\frac{1}{36}+\frac{\sqrt{3}}{3\pi}+\frac{3}{\pi^{2}}\approx 0.516.$ (10)

This is shown by the feint grey horizontal line in figure 3.

##  3 Magneto-frictional results

For a one-dimensional magnetic field of the form $B_{x}=B_{z}=0$,
$B_{y}=B_{y}(x,t)$, the frictional velocity (2) reduces to

$u_{x}=-\frac{1}{\nu}\frac{\upartial{}}{\upartial{x}}\left(\frac{B_{y}^{2}}{2\mu_{0}}\right),\hskip
28.45274ptu_{y}=u_{z}=0,$ (11e,f)

so equation (1) takes the form

$\frac{\upartial{B_{y}}}{\upartial{t}}=\frac{1}{\mu_{0}}\frac{\upartial{}}{\upartial{x}}\left(\frac{B_{y}^{2}}{\nu}\frac{\upartial{B_{y}}}{\upartial{x}}\right).$
(12)

The behaviour of this equation depends on the chosen form of the friction
coefficient $\nu$. Overall scaling of $\nu$ is simply equivalent to scaling
the time variable, so we only need to consider its functional form. In the
following subsections, we will consider three illustrative cases:

  1. 1.

$\nu=B_{y}^{2}$;

  2. 2.

$\nu=1$;

  3. 3.

$\nu=B_{y}^{2}+\varepsilon\mathrm{e}^{-B_{y}^{2}/\varepsilon}$, where
$\varepsilon$ is a small constant ($10^{-1}$ in our computations).

These cover the basic forms that have been used by different authors in the
literature. The results of these three test cases are shown in figure 4, and
discussed in the following subsections. For cases (ii) and (iii), equation
(12) was solved with the Crank-Nicolson method, treating the nonlinearity by
Picard iteration (see Appendix 5 for details). Results shown in figure 4 used
mesh resolution $n_{x}=1024$. The unsigned magnetic flux in figure 4(e) is
defined as

$|\varPhi|=\int_{-1}^{1}|B_{y}|\,\mathrm{d}x.$ (13)

Figure 4: MF results for the test case with each form of $\nu$. Panels (a) to
(d) show the solutions at selected times, according to the legend in panel
(f). The thick gray line shows the minimum energy solution from section 2.2.
Panel (e) shows the total unsigned magnetic flux as a function of time, and
panel (f) shows the magnetic energy. Vertical dashed lines in (e) and (f)
indicate the times of panels (a)-(d), and vertical solid lines in (e) indicate
the breakdown times $t^{*}$ and $t^{*}_{\varepsilon}$ as described in the
text. We take $\varepsilon=10^{-1}$. (Colour online)

###  3.1 Case (i) – linear diffusion

As noted in section 1, applications of MF have commonly set $\nu=B^{2}$, in
order to speed up relaxation in regions of low $B$. In our one-dimensional
case, this reduces (12) to the linear diffusion equation

$\frac{\upartial{B_{y}}}{\upartial{t}}=\frac{1}{\mu_{0}}\frac{\upartial{}^{2}B_{y}}{\upartial{x^{2}}},$
(14)

whose solution remains smooth for all time. For the test case (1), the
solution is precisely $B_{y}(x,t)=1/2+\sin(\pi x)\exp(-\pi^{2}t)$, shown by
the green dotted lines in figure 4.

Unfortunately, equation (14) no longer conserves (unsigned) magnetic flux, as
evidenced by the green dotted line in figure 4(e). (After approximately
$t=0.07$ unsigned flux is conserved because the two null points have been
completely eliminated.) Thus it violates a fundamental property of ideal
relaxation. Although we are still solving the induction equation (5d), there
is no contradiction because the frictional velocity from (3) is

$u_{x}=-\frac{1}{\mu_{0}B_{y}}\frac{\upartial{B_{y}}}{\upartial{x}},$ (15)

which is singular at null points so that flux conservation breaks down (e.g.
Wilmot-Smith et al., 2005). Not surprisingly, the magnetic energy of this
solution in figure 4(f) reaches far below that of the ideal minimum energy
state (in the limit, it tends to $1/4$, approximately half that of the ideal
minimum).

More generally, it would be only the perpendicular current that would be
diffused with $\nu=B^{2}$. To see this, note that, in three dimensions,
setting $\nu=B^{2}$ in equations (1) and (2) would give

$\frac{\upartial{\bm{B}}}{\upartial{t}}=-{\bm\nabla}\times\bm{J}_{\perp},\hskip
28.45274pt\bm{J}_{\perp}=\bm{J}-\frac{\bm{J}{\bm\cdot}\bm{B}}{B^{2}}\bm{B}.$
(16d,e)

Nevertheless, our one-dimensional test suggests that the method is unable in
general to accurately reconstruct an ideal relaxed state if null points are
present. For example, we propose that this diffusion, rather than numerical
error, could have been responsible for the change in magnetic topology
observed in two dimensions by Yang et al. (1986).

###  3.2 Case (ii) – ambipolar diffusion

Next, consider the case of constant $\nu$. In that case, equation (12) is
equivalent to the induction equation for a partially-ionised plasma under
ambipolar diffusion and vanishing plasma velocity (Brandenburg and Zweibel,
1994; Hoyos et al., 2010). It is a particular case of the porous medium
equation (Vazquez, 2006). An important property of this equation is that
discontinuous current sheets develop at finite time from null points $B_{y}=0$
in the initial state. We can predict the time of this breakdown by considering
the slope $w_{x_{0}}(t)=(\upartial B_{y}/\upartial x)_{x_{0}}$ at a null point
$x=x_{0}$ (Low, 2013). As long as the solution remains smooth, the null point
cannot move because $u_{x}(x_{0},t)=0$ by equation (3). So differentiating
(12) with respect to $x$ and setting $B_{y}=0$ shows that

$\frac{\mathrm{d}w_{x_{0}}}{\mathrm{d}t}=\frac{2}{\mu_{0}}w_{x_{0}}^{3}$ (17)

and hence the slope at the null, $w_{x_{0}}(t)$, obeys

$w_{x_{0}}^{2}(t)=\frac{w_{x_{0}}^{2}(0)}{1-4w_{x_{0}}^{2}(0)t}.$ (18)

It follows that the smooth solution will break down and a current sheet will
form at time $t^{*}=[4w_{x_{0}}^{2}(0)]^{-1}$. Similar behaviour is expected
in higher dimensions, and indeed finite-time formation of a discontinuous
current sheet at a three-dimensional null point is strongly suggested by the
numerical results of Pontin and Huang (2012), who also use MF with constant
$\nu$. Such breakdown of smoothness is not in itself a problem in our test
case, where we know that the minimum energy state under ideal relaxation
neglecting fluid pressure would have discontinuous current sheets.
Incidentally, Craig and Litvinenko (2005) show that adding a pressure gradient
term to the magneto-frictional velocity (2) can prevent the formation of
discontinous current sheets in one-dimensional configurations (see their
Appendix A). But in higher dimensions, Pontin and Craig (2005) demonstrate
that adding a pressure term is insufficient to prevent current-sheet collapse,
since the Lorentz force is not generally irrotational.

After $t^{*}$, which takes the value $t^{*}=(3\pi^{2})^{-1}$ for our test
case, we can continue to follow the evolution as a weak solution, since (12)
has the form of a conservation law. However, two properties that held while
the solution was smooth are no longer true: the null points can now move, and
the unsigned magnetic flux is no longer conserved. We now discuss each of
these in turn.

Figure 5: Verification that the numerical MF solution with $\nu=1$ satisfies
the formulae (19) and (20b). The solid curves show (a) the position of the
right-hand null/current sheet and (b) the total unsigned magnetic flux. The
dashed curves show predictions of these quantities found by estimating the
jumps in (19) and (20b) numerically and integrating these expressions in time.
(Colour online)

It is evident from the time snapshots in figure 4(a–d) that, once the null
points have degenerated into discontinuous current sheets, these start to
move. Their speed of movement is given by the standard Rankine-Hugoniot
condition,

$\frac{\mathrm{d}x_{0}}{\mathrm{d}t}=-[B_{y}]^{-1}_{x_{0}(t)}\left[\frac{B_{y}^{2}}{\nu}\frac{\upartial{B_{y}}}{\upartial{t}}\right]_{x_{0}(t)},$
(19)

where square brackets denote the limiting jump in a quantity across the
current sheet $x=x_{0}(t)$. Note that the right-hand side would vanish if the
ideal minimum energy state were reached, although it is not actually reached.
Figure 5(a) shows that our numerical solution recovers the correct speed for
the current sheet that starts at $x_{0}(t^{*})=-1/6$. It is clear that the
objection of Low (2013) to MF – that null points cannot move – holds only
while the solution remains smooth, and is not relevant in the long term.

Next we consider the breakdown in conservation of unsigned flux after $t^{*}$.
This is clearly evident in figure 4(e), and is not caused by numerical error.
The fact that ambipolar diffusion leads to loss of magnetic flux after the
formation of current sheets was shown by Hoyos et al. (2010) for a $B_{y}$
profile with a single step. They called it “reconnection in the absence of
ohmic diffusivity”. For our periodic domain with discontinuous current sheets,
the dissipation rate of unsigned flux may be expressed in terms of jumps at
the current sheets. In our example with two current sheets $x_{0}(t)$ and
$x_{0}^{\prime}(t)$, careful evaluation of the integral (13) shows – after
some algebra – that

$\displaystyle\frac{\mathrm{d}|\varPhi|}{\mathrm{d}t}$
$\displaystyle=[B_{y}]^{-1}_{x_{0}^{\prime}(t)}\left[\frac{2B_{y}^{3}}{\nu}\frac{\upartial{B_{y}}}{\upartial{x}}\right]_{x_{0}^{\prime}(t)}-[B_{y}]^{-1}_{x_{0}(t)}\left[\frac{2B_{y}^{3}}{\nu}\frac{\upartial{B_{y}}}{\upartial{x}}\right]_{x_{0}(t)}$
(20a)
$\displaystyle=-4[B_{y}]^{-1}_{x_{0}(t)}\left[\frac{B_{y}^{3}}{\nu}\frac{\upartial{B_{y}}}{\upartial{x}}\right]_{x_{0}(t)},$
(20b)

where in the second line we used the symmetry of the setup. Although we cannot
evaluate these quantities analytically, equation (20b) gives us a relation
that we can use to check that the dissipation observed in the numerical
solution is accurate. Indeed, figure 5(b) shows that our numerical method
obeys this relation to reasonable accuracy, despite the difficulty in
estimating the limiting quantities on either side of the current sheet from
the numerical solution. This arises since $B_{y}$ is not uniform on each side
of the current sheet but varies like $(x-x_{0})^{1/3}$ (Brandenburg and
Zweibel, 1994). In conclusion, it is clear that the dissipation of unsigned
flux after $t^{*}$ is a real phenomenon. This clearly prevents the method from
recovering the correct minimum-energy state, although it is evident from
figure 4(e) that the unwanted dissipation is much slower than for linear
diffusion.

###  3.3 Case (iii) – linear diffusion with limiting

Figure 6: Functional forms of the friction coefficient with limiting, for
$\varepsilon=10^{-1}$, compared to the linear diffusion case. (Colour online)

The final friction coefficient we consider is
$\nu=B_{y}^{2}+\varepsilon\mathrm{e}^{-B_{y}^{2}/\varepsilon}$. This
functional form is plotted in figure 6, and is seen to limit the minimum value
of $\nu$ to $\varepsilon$, while keeping $\nu\approx B_{y}^{2}$ away from null
points, as with linear diffusion. Such a limit is typically (albeit tacitly)
employed in numerical calculations; for example, Yeates (2014) use the simpler
non-smooth form $\nu=\max\\{B_{y}^{2},\varepsilon\\}$, also shown in figure 6.
The parameter $\varepsilon$ would be set to a small number. The resulting
equation is

$\frac{\upartial{B_{y}}}{\upartial{t}}=\frac{1}{\mu_{0}}\frac{\upartial{}}{\upartial{x}}\left(\frac{B_{y}^{2}}{B_{y}^{2}+\varepsilon\mathrm{e}^{-B_{y}^{2}/\varepsilon}}\frac{\upartial{B_{y}}}{\upartial{x}}\right).$
(21)

This limited form of $\nu$ is interesting mathematically because it prevents
the frictional velocity (3) from being singular at null points, thus
preserving flux conservation. However, there is a catch: near null points, the
equation looks like the ambipolar diffusion case (constant $\nu$). Thus we
might expect the finite-time breakdown of smoothness and formation of
discontinuous current sheets at nulls. Indeed we can apply the argument of Low
(2013) to equation (21) to see – after some algebra – that a smooth solution
with a null point will break down at
$t^{*}_{\varepsilon}=\varepsilon[4w_{x_{0}}^{2}(0)]^{-1}$, where
$w_{x_{0}}(0)$ is again the initial slope of $B_{y}$ at the null. In
particular, for $\varepsilon<1$, this will occur more rapidly than for the
ambipolar diffusion case.

The orange dashed curves in figure 4 show the numerical solution with the
limiting form of $\nu$ for $\varepsilon=10^{-1}$, and indeed a discontinuous
current sheet forms already at $t=t^{*}_{\varepsilon}$. As for the constant
$\nu$ solution, the null points remain stationary and flux is conserved until
this breakdown, but both properties are violated thereafter. It should be
noted that a large value of $\varepsilon=10^{-1}$ was used here for
illustration. As $\varepsilon$ is reduced, $t^{*}_{\varepsilon}$ gets earlier
and earlier, while the current sheets at any fixed later time have less and
less of a jump in $B_{y}$, with the solution tending to the linear diffusion
case as $\varepsilon\to 0$. For the small $\varepsilon$ typically used in
simulations to avoid division by zero, the solution would be essentially
indistinguishable from the linear diffusion case, despite the fact that the
frictional velocity is initially non-singular.

##  4 Discussion

We have shown that, when a magnetic null point is present, magneto-frictional
(MF) relaxation is unable to accurately predict the relaxed state that would
be obtained by an ideal MHD evolution. We have demonstrated two separate
reasons for this. Firstly, as shown in section 2.1, a relaxation under the
full MHD equations starting from a uniform pressure will – at least in our
one-dimensional example – build up a non-trivial distribution of fluid
pressure, owing to the fact that magnetic energy is transformed to internal
energy. In particular, fluid pressure accumulates at the magnetic null points
in order to reach total pressure balance. Thus the true equilibrium will be
magnetostatic rather than force-free. In effect, the plasma-beta is high near
the null points, even if it is low elsewhere. This motivates the use of
relaxation schemes that include fluid pressure, while still avoiding the need
for small timesteps to resolve MHD waves. As discussed, one approach is to
modify the MF velocity to $\bm{u}=\nu^{-1}(\bm{J}\times\bm{B}-{\bm\nabla}p)$.
To determine $p$ in two-dimensional simulations, Linardatos (1993) used the
requirement that ${\bm\nabla}{\bm\cdot}\bm{u}=0$ (see also Moffatt and Dormy,
2019). Unfortunately, this would not work for our one-dimensional test case,
where the relaxation velocity is necessarily compressible. Instead, one could
write ${\bm\nabla}p=\beta{\bm\nabla}\rho$, and evolve $\rho$ by mass
conservation (Craig and Litvinenko, 2005; Candelaresi et al., 2015).

A more significant limitation of MF without pressure is that it fails to
respect conservation of unsigned magnetic flux, which would follow from the
ideal induction equation (5d) provided that $\bm{u}$ remains smooth. The
problem is that $\bm{u}$ does not remain smooth for any of the forms of MF in
common use. When $\nu=B_{y}^{2}$, our one-dimensional problem reduces to a
linear diffusion equation. When $\nu$ is constant, or becomes constant near to
a null point, then $\bm{u}$ is initially smooth but this breaks down at a
finite time. In itself, we would argue (in contrast to Low, 2013) that this
formation of discontinuous current sheets is not a problem, since the state of
minimum magnetic energy when fluid pressure is neglected is a discontinuous
one (section 2.2). However, the ensuing weak solution obtained does not
conserve unsigned flux, so does not evolve toward this expected solution. For
one-dimensional configurations such as that investigated here, Craig and
Litvinenko (2005) show that adding plasma pressure can also prevent the
formation of a discontinuous current sheet, and hence also the breakdown of
flux conservation.

Figure 7: Result of viscous relaxation (with $\mu=1$) for the test case. As in
figure 4, panels (a) and (b) show snapshots of the solution during relaxation,
while (c) and (d) show the unsigned magnetic flux and magnetic energy,
respectively. The thick gray line shows the minimum energy solution from
section 2.2. (Colour online)

To end this paper on a more positive note, we point out an alternative method
that is able to avoid the second limitation (flux non-conservation) and reach
the expected minimum energy state from section 2.2. This is the viscous
relaxation scheme (Bajer and Moffatt, 2013) where we again solve the induction
equation (1), but determine $\bm{u}$ from the force balance

$\mu_{s}\Delta\bm{u}+\left(\frac{1}{3}\mu_{s}+\mu_{b}\right){\bm\nabla}({\bm\nabla}{\bm\cdot}\bm{u})+\bm{J}\times\bm{B}=\bm{0}.$
(22)

In our one-dimensional case this reduces to

$\mu\frac{\upartial{}^{2}u_{x}}{\upartial{x}^{2}}-\frac{\upartial{}}{\upartial{x}}\left(\frac{B_{y}^{2}}{2\mu_{0}}\right)=0,$
(23)

where $\mu=4\mu_{s}/3+\mu_{b}$. For our periodic system, this scheme again has
the property that the magnetic energy decreases monotonically, with

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{V}\frac{B^{2}}{2\mu_{0}}\,\mathrm{d}V=-\mu\int_{V}\left|\frac{\upartial{u_{x}}}{\upartial{x}}\right|^{2}\,\mathrm{d}V.$
(24)

In the periodic system, equation (23) specifies $u_{x}$ only up to an additive
constant. In view of the symmetry of our test problem, we set
$u_{x}(-1/2,t)=0$, so solving (23) for $u_{x}$ gives

$u_{x}(x,t)=\frac{1}{\mu}\int_{-1/2}^{x}\frac{B_{y}^{2}}{2\mu_{0}}\,\mathrm{d}x-\frac{x+1/2}{2\mu}\int_{-1}^{1}\frac{B_{y}^{2}}{2\mu_{0}}\,\mathrm{d}x.$
(25)

Substituting this into the induction equation (5d) gives an integro-
differential evolution equation for $B_{y}$. Figure 7 shows a numerical
solution on the same uniform mesh used for the MF simulations, but using a
simple upwind scheme (as the equation is now hyperbolic). Unlike MF, we see
that the viscous relaxation method evolves rapidly toward the minimum energy
state. In future, it would be fruitful to explore as an alternative to MF in
the context of the solar coronal magnetic field, where null points are
prevalent.

## Acknowledgements

This work was supported by UKRI/STFC research grant ST/S000321/1. The author
thanks the two anonymous referees for improving the paper.

## References

  * Bajer and Moffatt (2013) Bajer, K. and Moffatt, H.K., Magnetic relaxation, current sheets, and structure formation in an extremely tenuous fluid medium. Astrophys. J., 2013, 779, 169(14pp).
  * Bhowmik and Yeates (2021) Bhowmik, P. and Yeates, A.R., Two classes of eruptive events during solar minimum. Solar Phys., 2021, 296, 109(30pp).
  * Brandenburg and Zweibel (1994) Brandenburg, A. and Zweibel, E.G., The formation of sharp structures by ambipolar diffusion. Astrophys. J. Lett., 1994, 427, L91–L94.
  * Candelaresi et al. (2015) Candelaresi, S., Pontin, D.I. and Hornig, G., Magnetic field relaxation and current sheets in an ideal plasma. Astrophys. J., 2015, 808, 134(8pp).
  * Cheung and DeRosa (2012) Cheung, M.C.M. and DeRosa, M.L., A method for data-driven simulations of evolving solar active regions. Astrophys. J., 2012, 757, 147(10pp).
  * Chodura and Schlüter (1981) Chodura, R. and Schlüter, A., A 3D code for MHD equilibrium and stability. J. Comp. Phys., 1981, 41, 68–88.
  * Craig and Litvinenko (2005) Craig, I.J.D. and Litvinenko, Y.E., Current singularities in planar magnetic X points of finite compressibility. Phys. Plasmas, 2005, 12, 032301(10pp).
  * Craig and Sneyd (1986) Craig, I.J.D. and Sneyd, A.D., A dynamic relaxation technique for determining the structure and stability of coronal magnetic fields. Astrophys. J., 1986, 311, 451–459.
  * Goldstraw et al. (2018) Goldstraw, E.E., Hood, A.W., Browning, P.K. and Cargill, P.J., Comparison of methods for modelling coronal magnetic fields. Astron. Astrophys., 2018, 610, A48(15pp).
  * Guo et al. (2016) Guo, Y., Xia, C., Keppens, R. and Valori, G., Magneto-frictional modeling of coronal nonlinear force-free fields. I. Testing with analytic solutions. Astrophys. J., 2016, 828, 82(21pp).
  * Hesse and Birn (1993) Hesse, M. and Birn, J., Three-dimensional magnetotail equilibria by numerical relaxation techniques. J. Geophys. Res., 1993, 98, 3973–3982.
  * Hoeksema et al. (2020) Hoeksema, J.T., Abbett, W.P., Bercik, D.J., Cheung, M.C.M., DeRosa, M.L., Fisher, G.H., Hayashi, K., Kazachenko, M.D., Liu, Y., Lumme, E., Lynch, B.J., Sun, X. and Welsch, B.T., The Coronal Global Evolutionary Model: using HMI vector magnetogram and doppler data to determine coronal magnetic field evolution. Astrophys. J. Supp. Ser., 2020, 250, 28(15pp).
  * Hoyos et al. (2010) Hoyos, J.H., Reisenegger, A. and Valdivia, J.A., Asymptotic, non-linear solutions for ambipolar diffusion in one dimension. Mon. Not. R. Astron. Soc., 2010, 408, 1730–1741.
  * Klimchuk and Sturrock (1992) Klimchuk, J.A. and Sturrock, P.A., Three-dimensional force-free magnetic fields and flare energy buildup. Astrophys. J., 1992, 385, 344–353.
  * Linardatos (1993) Linardatos, D., Determination of two-dimensional magnetostatic equilibria and analogous Euler flows. J. Fluid Mech., 1993, 246, 569–591.
  * Low (2013) Low, B.C., Newtonian and non-newtonian magnetic-field relaxations in solar-coronal MHD. Astrophys. J., 2013, 768, 7(17pp).
  * Lowder and Yeates (2017) Lowder, C. and Yeates, A., Magnetic flux rope identification and characterization from observationally driven solar coronal models. Astrophys. J., 2017, 846, 106(13pp).
  * Mackay et al. (2011) Mackay, D.H., Green, L.M. and van Ballegooijen, A., Modeling the dispersal of an active region: quantifying energy input into the corona. Astrophys. J., 2011, 729, 97(11pp).
  * Mackay and Yeates (2012) Mackay, D.H. and Yeates, A.R., The sun’s global photospheric and coronal magnetic fields: observations and models. Living Rev. Solar Phys., 2012, 9, 6(63pp).
  * Meyer and Mackay (2016) Meyer, K.A. and Mackay, D.H., Modeling the sun’s small-scale global photospheric magnetic field. Astrophys. J., 2016, 830, 160(13pp).
  * Moffatt and Dormy (2019) Moffatt, H.K. and Dormy, E., Self-exciting Fluid Dynamos, 2019 (Cambridge University Press).
  * Pomoell et al. (2019) Pomoell, J., Lumme, E. and Kilpua, E., Time-dependent data-driven modeling of active region evolution using energy-optimized photospheric electric fields. Solar Phys., 2019, 294, 41(28pp).
  * Pontin and Craig (2005) Pontin, D.I. and Craig, I.J.D., Current singularities at finitely compressible three-dimensional magnetic null points. Phys. Plasmas, 2005, 12, 072112(8pp).
  * Pontin and Huang (2012) Pontin, D.I. and Huang, Y.M., On the formation of current sheets in response to the compression or expansion of a potential magnetic fField. Astrophys. J., 2012, 756, 7(6pp).
  * Roumeliotis (1996) Roumeliotis, G., The “stress-and-relax” method for reconstructing the coronal magnetic field from vector magnetograph data. Astrophys. J., 1996, 473, 1095–1103.
  * Stone et al. (2008) Stone, J.M., Gardiner, T.A., Teuben, P., Hawley, J.F. and Simon, J.B., Athena: a new code for astrophysical MHD. Astrophys. J. Supp. Ser., 2008, 178, 137–177.
  * Toriumi et al. (2020) Toriumi, S., Takasao, S., Cheung, M.C.M., Jiang, C., Guo, Y., Hayashi, K. and Inoue, S., Comparative study of data-driven solar coronal field models using a flux emergence simulation as a ground-truth data set. Astrophys. J., 2020, 890, 103(13pp).
  * Valori et al. (2005) Valori, G., Kliem, B. and Keppens, R., Extrapolation of a nonlinear force-free field containing a highly twisted magnetic loop. Astron. Astrophys., 2005, 433, 335–347.
  * Valori et al. (2010) Valori, G., Kliem, B., Török, T. and Titov, V.S., Testing magnetofrictional extrapolation with the Titov-Démoulin model of solar active regions. Astron. Astrophys., 2010, 519, A44(14pp).
  * van Ballegooijen et al. (2000) van Ballegooijen, A.A., Priest, E.R. and Mackay, D.H., Mean field model for the formation of filament channels on the sun. Astrophys. J., 2000, 539, 983–994.
  * Vazquez (2006) Vazquez, J.L., The Porous Medium Equation, 2006 (Clarendon Press).
  * Wiegelmann and Sakurai (2021) Wiegelmann, T. and Sakurai, T., Solar force-free magnetic fields. Living Rev. Solar Phys., 2021, 18, 1(67pp).
  * Wilmot-Smith et al. (2005) Wilmot-Smith, A.L., Priest, E.R. and Hornig, G., Magnetic diffusion and the motion of field lines. Geophys. Astrophys. Fluid Dyn., 2005, 99, 177–197.
  * Yang et al. (1986) Yang, W.H., Sturrock, P.A. and Antiochos, S.K., Force-free magnetic fields: the magneto-frictional method. Astrophys. J., 1986, 309, 383–391.
  * Yardley et al. (2021) Yardley, S.L., Mackay, D.H. and Green, L.M., Simulating the coronal evolution of bipolar active regions to investigate the formation of flux ropes. Solar Phys., 2021, 296, 10(27pp).
  * Yeates (2014) Yeates, A.R., Coronal magnetic field evolution from 1996 to 2012: continuous non-potential simulations. Solar Phys., 2014, 289, 631–648.

##  5 Numerical method for MF

The MF equation (12) has the form of a nonlinear diffusion equation, and is
either parabolic or degenerate parabolic, depending on the form of $\nu$. As
such it is best treated with an implicit method, although solar physics
applications have often used an explicit method with the resulting timestep
restriction. In this paper we adopt an implicit Crank-Nicolson scheme,
treating the nonlinearity with Picard iteration. For shorthand, let $B=B_{y}$
and $\beta(B^{2})=B^{2}/[\mu_{0}\nu(B^{2})]$, so that the equation is

$\frac{\upartial{B}}{\upartial{t}}=\frac{\upartial{}}{\upartial{x}}\left(\beta(B^{2})\frac{\upartial{B}}{\upartial{x}}\right).$
(26)

We take a uniform mesh $x_{j}=j\Delta x-1$ for $j=0,\cdots,n_{x}-1$, with
spacing $\Delta x=2/(n_{x}-1)$. The magnetic field values
$B_{j+\frac{1}{2}}^{n}$ are located at cell centres
$\frac{1}{2}(x_{j}+x_{j+1})$, where $n$ denotes the time level $t_{n}=n\Delta
t$. A single timestep, to compute $B_{j+\frac{1}{2}}^{n+1}$ from
$B_{j+\frac{1}{2}}^{n}$, involves iteratively solving the Crank-Nicolson
formula over $k$, where

$\displaystyle-\alpha_{j}^{n,k}B_{j-\frac{1}{2}}^{n,k+1}+\left(1+\alpha_{j}^{n,k}+\alpha_{j+1}^{n,k}\right)B_{j+\frac{1}{2}}^{n,k+1}-\alpha_{j+1}^{n,k}B_{j+\frac{3}{2}}^{n,k+1}$
$\displaystyle\hskip
113.81102pt=\alpha_{j}^{n,0}B_{j-\frac{1}{2}}^{n}+\left(1-\alpha_{j}^{n,0}-\alpha_{j+1}^{n,0}\right)B_{j+\frac{1}{2}}^{n}+\alpha_{j+1}^{n,0}B_{j+\frac{3}{2}}^{n},\hskip
28.45274pt$ (27)

starting from $B_{j+\frac{1}{2}}^{n,0}=B_{j+\frac{1}{2}}^{n}$. The
coefficients comprise the scaled “diffusivity” $\beta$ evaluated at the mesh
points,

$\alpha_{j}^{n,k}=\frac{\Delta t}{2\Delta
x^{2}}\beta\left(\bigl{[}B^{2}\bigr{]}_{j}^{n,k}\right),$ (28a) in which
$\bigl{[}B^{2}\bigr{]}_{j}^{n,k}=\frac{1}{4}\biggl{[}\Bigl{(}B_{j-\frac{1}{2}}^{n}\Bigr{)}^{\\!2}+\Bigl{(}B_{j+\frac{1}{2}}^{n}\Bigr{)}^{\\!2}+\Bigl{(}B_{j-\frac{1}{2}}^{n,k}\Bigr{)}^{\\!2}+\Bigl{(}B_{j+\frac{1}{2}}^{n,k}\Bigr{)}^{\\!2}\biggr{]}.$
(28b)

Equation (27) is iteratively solved until
$\left\|B_{j+\frac{1}{2}}^{n,k+1}-B_{j+\frac{1}{2}}^{n,k}\right\|_{\infty}<10^{-14}$,
then we set $B_{j+\frac{1}{2}}^{n+1}=B_{j+\frac{1}{2}}^{n,k+1}$.

We have verified the convergence against exact solutions for (i) linear
diffusion when $\nu=B^{2}$, and (ii) ambipolar diffusion when $\nu=1$. For
(ii) we do not have the exact solution for our test problem, so we tested
instead with the well-known Barenblatt solution (Vazquez, 2006) in the form

$B_{y}(x,t)=\sqrt{\max\left\\{0,\frac{1}{10}\tilde{t}^{-1/2}-\frac{x^{2}}{12}\tilde{t}^{-1}\right\\}},\hskip
28.45274pt\tilde{t}=\dfrac{1}{3}\left(t+\dfrac{1}{10}\right).$ (29)

For case (i), and in smooth regions of case (ii), the expected second-order
convergence is recovered, where we took $\Delta t=10\Delta x^{2}$. The
position of the moving “front” at $x=(\sqrt{6/5})\tilde{t}^{1/4}$, where the
slope of $B_{y}$ in the exact solution becomes infinite and discontinuous,
still converges, but only to first order.