Śniegowska et al

*Corresponding author.

# Modeling changing-look (CL) AGN phenomenon in 1D using accretion disk
instabilities

Marzena Śniegowska  Mikołaj Grzędzielski  Bożena Czerny  Agnieszka Janiuk
Nicolaus Copernicus Astronomical Center (PAN), ul. Bartycka 18, 00-716 Warsaw,
Poland  Center for Theoretical Physics, Polish Academy of Sciences, Al.
Lotników 32/46, 02-668 Warsaw, Poland<EMAIL_ADDRESS>Śniegowska M
Grzędzielski M  Czerny B  Janiuk A

###### Abstract

Apart from regular, low-level stochastic variability, some AGN occasionally
show exceptionally large changes in the luminosity, spectral shape, and/or
X-ray absorption. The most notable are the changes of the spectral type when
the source classified as a Seyfert 1 becomes a Seyfert 2 galaxy or vice versa.
Thus a name was coined of ’Changing-Look AGN’ (CL AGN). The origin of this
phenomenon is still unknown, but for most of the sources, there are strong
arguments in favor of the intrinsic changes.

Understanding the nature of such rapid changes is a challenge to the models of
black hole accretion flows since the timescales of the changes are much
shorter than the standard disk viscous timescales. We aim to model the CL AGN
phenomenon assuming that the underlying mechanism is the time-dependent
evolution of a black hole accretion disk unstable due to the dominant
radiation pressure. We use a 1-dimensional, vertically integrated disk model,
but we allow for the presence of the hot coronal layer above the disk and the
presence of the inner purely hot flow. We focus on the variability timescales
and amplitudes, which can be regulated by the action of large-scale magnetic
fields, the description of the disk-corona coupling and the presence of an
inner optically thin flow, like Advection-Dominated Accretion Flow (ADAF). We
compare model predictions for the accretion disk around black hole mass 107M⊙.

###### keywords:

accretion, accretion disks, viscosity

††articletype: Article Type

, , , and (2021), Modelling changing-look (CL) AGN phenomenon in 1D using
accretion disk instabilities, Astronomische Nachrichten, ???.

##  1 INTRODUCTION

In the simplest theoretical framework, active galactic nuclei (AGN) are
sources in which a supermassive black hole is surrounded by a geometrically
thin, optically thick accretion disk [Shakura Sunyaev (1973)]. For some
objects, though, modifications to accretion disk, like advection-dominated
accretion flow - ADAF or Slim Accretion Disks (for review see [Czerny
(2019)]), are needed to better represent the observed spectra. AGN also show
variability in different timescale, pattern and wavelengths like X-ray
(Lawrence Papadakis, 1993), optical (Sesar ., 2007) and infrared (Kozłowski .,
2016) bands.

Generally, in accretion disk heating processes is mediated by the turbulent
magnetic field. Thus the stochastic variability observed in light curves of
AGN is not surprising. However, objects with significant changes in the
luminosity and the noticeable (dis)appearance of broad emission lines are less
likely to be explained by stochastic variability. They are called changing-
look (CL) AGN and the discussion about their physical nature is coming under
scrutiny in recent years.

In this paper, we focus on accretion disk around supermassive black holes
(107M⊙). The variations in global accretion rate in this kind of accretion
disks are one of the candidates for CL AGN physical explanation. However, it
seems that those objects may be explained by different scenarios and just one
explanation for them is not likely. CL AGN differ between each other in the
shape of the light curve, timescale of changes, and its amplitude. For
example, changes in some objects (like Mrk 590 Kokubo Minezaki (2020)) can be
explained by the extrinsic matter which causes the dimming of the flux and
spectral changes. Whereas, for others (like in sources investigated by
(Hutsemékers ., 2019)) the obscuration is ruled out because of no changes in
polarisation degree, which should appear if we assume dusty matter is
temporary obscuring the source. Thus, intrinsic changes are in favor to
explain it in most cases. Intrinsic changes in accretion disks are rather a
broad concept and may include scenarios with tidal disruption events (i.e SDSS
J224113-012108 investigated by Zhang (2021)), tidal interaction between disks
in supermassive black hole binaries (Wang Bon, 2020), the presence and
disappearance of the warm corona (Noda Done, 2018), changes in the
magnetization of accretion disk Scepi . (2021) or magnetic accretion disk-
outflows (Feng ., 2021). The rapidly increasing number of observations of CL
AGN (Sánchez-Sáez ., 2021; Graham ., 2020; MacLeod ., 2016) and variety of
possible scenarios of changes motivates us to explore this field.

To model CL AGN phenomenon (especially for objects with repetitive outbursts),
we assume the radiation pressure instability mechanism. In the current work,
we use time-evolution code GLADIS (Global Accretion Disk Instability
Simulation) developed by Janiuk . (2002) (recently publicly available version
is presented in Janiuk (2020)), which allows performing a time-dependent
simulation of accreting matter onto massive objects. This code was already
used to model the ’heartbeat’ outbursts in microquasars (e.g. for IGR
J17091-3624 by Janiuk ., 2015) as caused by radiation pressure instabilities.
The idea of explaining the changes in CL AGN with the radiation pressure
instability was used in Sniegowska . (2020), however, the concept of the
geometry of the system was very simplistic, reducing the unstable region to a
single zone both in radial and vertical direction. In this work, we explore
the instability using 1-D radially resolved model. In Fig. 1 we present the
scheme of the model which we use. This model contains the Advection-Dominated
Accretion Flow (ADAF) in which hot coronal matter is freely advected inwards,
the unstable standard cold disk zone radially resolved, the stable outer disk
and the geometrically thick, optically thin corona above the disk. Therefore,
in comparison to Sniegowska . (2020), this model is extended by two important
aspects: (i) the presence of the disk corona, and (ii) radially resolved
unstable cold disk. The grid in the unstable zone is improving the model,
because the ring in which instabilities are present is not artificially
limited in its size. Disk instabilities in system disk-corona interacting were
investigated for the microquasar GRS 1915+105 by Janiuk Czerny (2005), with
the disk/corona mass exchange model developed later on Janiuk Czerny (2007)
for variability in X-ray binaries in their soft state and AGN. The presence of
the inner ADAF was not included in these papers. The structure of the paper is
as follows: In Section 2 we introduce the physical assumptions which we use
throughout this work for 3 computed models. In Section 3 we show and describe
the preliminary cycle-behavior time-dependent model results. Finally, in
Section 4, we conclude on the current results and comment on the plans for
future work.

##  2 THEORY

In the standard accretion disk theory, (Shakura Sunyaev, 1973) the assumption
of the radiation pressure dominance leads to the accretion disk instability
(thermal and viscous) which was shown by Pringle . (1973) and Lightman Eardley
(1974). Those instabilities may appear in the broad range of black hole
masses’ objects.

Outbursts that may appear due to radiation pressure instabilities in
microquasars (masses $\approx$ 10 M⊙) were explored by Janiuk . (2002). This
kind of instabilities in objects with central black hole mass up to 108 M⊙ was
further described by Grzȩdzielski . (2017). The latter work gives also the
quantitative prescriptions for estimating periods and amplitudes of outbursts.
We use those prescriptions to discuss the predicted limits of the minimum
period for black hole mass of 107 M⊙.

Figure 1: The schematic view of the model sliced through the disk plane: from
the left black hole marked as a black circle, and inner hot ADAF, disk part,
which is unstable due to radiation pressure instability marked in blue, with
black stripes to emphasize radial grid, and outer cold stable disk represented
as a yellow area. The corona which is covering the disk is shown in magenta.

###  2.1 S-curve as the signpost for the initial conditions

A stability curve, in general, is a sequence of stationary solutions to
equations of disk structure equations, solved locally for surface density and
temperature. Before we perform the time-dependent simulation with the GLADIS
code, we compute the stability S-curves in an effective temperature - surface
density plane to get an overall picture (see Fig. 2). From the shape of the
S-curve, and specifically the range of its unstable branch (i.e., the negative
slope of $T-\Sigma$ dependence) we infer the information about the amplitudes
and timescales of outbursts. Along the unstable branch, the higher accretion
rates correspond to smaller surface densities. The position of the curve
depends on the radius. Hence, with changing the distance from the central
source, at a fixed temperature we switch between the stable and unstable
branch. Note, however, that the local effective temperature scales with the
radius and the local accretion rate, hence at a fixed accretion rate the
(sufficiently large) outer disk is always located at a stable branch. By
investigating S-curves behavior before running the time evolution code we
predict proper initial conditions and model parameters which allow us to
obtain instability of a desired range.

###  2.2 Set-up of the Model

The two time-dependent equations adopted by GLADIS code are equations of
hydrodynamics for surface density and temperature evolution under assumptions
of energy conservation, viscous diffusion, radiation pressure dominance (see
eq. 3.1 and 3.2 in Janiuk (2020))

In this work, we investigate 3 possible scenarios of mass exchange and
evaporation

  * •

mass exchange between disk and corona is on, and disk evaporation is due to
the flux generated in the disk and the corona - case A (see Janiuk Czerny
2005, their Eq. 17, with coefficients $B_{1}=0.5$, $B_{2}=0.5$). This
description represents the disk evaporation due to the electron conduction
between the disk and the corona,

  * •

mass exchange between disk and corona is on, and disk evaporation is due to
magnetic field - case B (see Janiuk Czerny 2007, their Eq. 10, their case B).
This description represents the corona heated by the magnetic buoyancy,

  * •

mass exchange between disk and corona is off, which implies no disk
evaporation - case C.

For all calculated models in this work, we assume that the heating term is
proportional to the total pressure, as was assumed by Shakura Sunyaev (1973).

We fix the following physical parameters throughout all the computed models:
the black hole mass 107M⊙, viscosity coefficient $\alpha$ = 0.02 for disk and
corona after Grzȩdzielski . (2017), global accretion rate as 0.007 solar
masses per year, the inner radius of the disk 20 RSchw, the outer radius of
the disk 100 RSchw00footnotetext: due to limitations in computing time we
decided to keep outer radius in this position, however, unstable zones are
more extended (see Janiuk Czerny (2011)). We discuss this issue later. Corona
temperature is fixed at the virial temperature, whereas the temperature of the
disc evolves in time. The total duration of each simulation is 274 years in
the source time.

##  3 RESULTS

In this section, we compare the properties of time-dependent models for
different underlying assumptions. We calculate 3 solutions based on different
assumptions about the disk/corona mass exchange (see Sec. 2.2).

Figure 2: Local stability curves for black hole mass 10${}^{7}M_{\odot}$,
$\alpha$ = 0.02 and different radii from the top: 5$R_{Schw}$ marked as red
curve, 25$R_{Schw}$ marked as blue curve and 100$R_{Schw}$ marked as green
curve.

###  3.1 The outburst cycle

The example of outburst for the object with the black hole mass 107M⊙ is shown
in the Fig 3. The luminosities (for disk and for corona) presented in this
work are the integrated radiation flux between inner and outer radius.

Each outburst has three characteristic phases:

  1. 1.

viscous phase (red part of the lightcurve), in which the approximate thermal
equilibrium holds, and the surface density increases till the system meets the
loss of equilibrium (lower branch in the S-curve),

  2. 2.

heating phase and hot viscous phase (the grey part of the lightcurve), when
the temperature in the inner regions of the disk increases rapidly, and the
accretion rate is enhanced, the system then evolves along the upper S-curve
branch till the turning point of the S-curve is met,

  3. 3.

cooling phase (blue part of the lightcurve), where the temperature and
accretion rate drop-down till the lower branch of the S-curve is achieved.

Outbursts can be then characterized by the amplitude, the period, and the
relative duration of the outburst in comparison with the total period, as
discussed by Grzȩdzielski . (2017) for their set of models.

Figure 3: The exemplary fragment of the disk lightcurve with cycle-steps
marked for black hole mass 10${}^{7}M_{\odot}$, Rin = 20Rshw and Rout =
100Rshw for case B. The grey part of the light curve represents heating phase,
the blue one advective phase and the red one diffusive phase. Note that we
present in this plot just part of the diffusive phase, not the full cycle.

Figure 4: Comparison of the modeled variations, for disc (gray line) and
corona (red dotted line) luminosity time evolution for 107M$\odot$ for three
scenarios: case A - mass exchange between disk and corona is on and mass
evaporation is due to the electron conduction (left panel), case B - mass
exchange between disk and corona is on and mass evaporation is due to magnetic
buoyancy (middle panel), case C - mass exchange between disk and corona is off
(right panel).

Figure 5: The exemplary fragment of the disk (marked in gray) and corona
(marked as dotted red line) lightcurve for case B.

###  3.2 Time-dependent solutions

In Fig.4 we show the modeled lightcurves of the disk (marked in gray) and the
corona (marked as red dotted line). This calculation is representative for
objects with observed repetitive outbursts (like NGC 1566, with log(M) = 6.92
$M_{\odot}$ (Woo Urry, 2002)). We fixed the viscosity coefficient in the disk
and corona $\alpha$ = 0.02 for both, as we stress out in Section 2.2. We
concentrate on the comparison of the three physical scenarios included in the
model. The lightcurve for corona is shown only in cases when disk-corona mass
exchange was used in computations. In the first model (left panel of Fig. 4 \-
case A), for assumed physical parameters disk and corona are stabilized. No
outbursts are seen, and the system simply evolves from the initial conditions
to the equilibrium solution. For two other investigated scenarios outbursts
are present. The characteristic three phases of the outbursts are thus visible
both when the disk/corona coupling is described by magnetic buoyancy (middle
panel of Fig. 4 \- case B) as well as for model with no corona included (right
panel of Fig. 4 \- case C). The diffusive phase is the longest of the three
phases.

The outbursts in both cases are very strong, with the amplitude (the ratio of
the maximum to the minimum flux in the disk $\frac{L_{max}}{L_{min}}$)
$\approx 10^{4}$. The period of the outburst for the case B is 156 years.
Flare duration to period ratio 00footnotetext: As a flare duration we assume
the time between $(L_{max}+L_{min})/2$ approached during the rise of the
outburst, and $(L_{max}+L_{min})/2$ approached during the dimming, following
parametrization from Grzȩdzielski . (2017). is 0.006.

For corona in case B (see Fig. 5 for the expanded illustration of both disk
and corona outburst) the luminosity starts to increase $\approx 0.5$ year
before the disk’s outburst and gets to the lowest luminosity level $\approx
0.5$ year before the disk’s dimming, so the period of the outburst is the same
for corona and accretion disk. The relative amplitude, in this case, is 50,
and the flare duration to period ratio is 0.009. Thus for these specific model
parameters the outbursts should be seen predominantly in the UV band, with
more moderate outbursts in X-ray. In case C the relative amplitude of the
outburst is 5000, by a factor 2 smaller than for case B, with the period of
147 years, which is shorter than in case B. Flare duration to period ratio is
0.007.

##  4 DISCUSSION

We performed simulations of the disk/corona/inner ADAF system (see Fig. 1)
under the radiation pressure instability in order to check if this mechanism
can explain the CL behavior, as proposed by Sniegowska . (2020). In the
preliminary one-zone model of Sniegowska . (2020), the expected timescales
were much shorter than the local viscous timescale due to the narrowness of
the instability zone. In the present computations, the outer radius is set at
100 $R_{Schw}$, and the extension of the instability zone is determined
numerically. During the outburst, the larger part of the disk becomes a
subject of instability than just the narrow unstable zone as described by the
stationary S-curves (see e.g. Grzȩdzielski . 2017). Therefore, the timescales
cannot be arbitrarily shortened by selecting the specific accretion rate and
the transition radius between the disk/corona part and the inner ADAF. Here we
show only one example for the global parameter values (black hole mass
$10^{7}M_{\odot}$, $\dot{m}=0.007$) but the same holds for other values of the
black hole and accretion rates (Sniegowska et al., in preparation).

We show that the model is actually sensitive to various assumptions about the
role of the electron conduction and/or magnetic field inside the disk, and the
description of the mass exchange between the disk and the corona. The disk is
stabilized in the model A. The outbursts are slightly shorter if the mass
exchange between the disk and the corona is turned off (see Fig. 4). The
duration of the outburst in the presented model is very short in comparison
with the separation of outbursts.

The timescales derived in Sect. 3.2 are the order of 150 years, shorter than
the outbursts timescales obtained by Grzȩdzielski . (2017) in case of modified
viscosity (using the maximum to the minimum luminosity in the outburst from
Fig. 4 and using with Eq. 31 of Grzȩdzielski . 2017 we obtain 7000 years).
This effect is not related to the presence of the corona which also transports
a fraction of the accreting material, but most of the effect is due to the
small outer disk radius used in current computations. Radiation pressure
instabilities mechanism would lead to outbursts of the much larger part of the
disk than 100$R_{Schw}$ (which we consider in this work as the outer radius of
simulated accretion disk). We checked that by computing one more model for
case B, adopting the outer radius at 80$R_{Schw}$, and those computations gave
a total period of 85 years, and apparently, the outer radius is the most
critical parameter of the current model.

This is a very interesting aspect of the model, but we currently cannot
address properly it without considerable modification of the GLADIS code.
Currently, we assume that the outer disk is stationary. If the disk outer
radius is indeed very small since the event is related to TDE, we would need
to model the whole disk as initially a ring, with a freely expanding outer
radius due to the viscous evolution. We will address this in the future.

Those results are not yet ready for detailed comparison with the lightcurve
for observed AGN, because the calculated zone is not extended enough to do so.
In the incoming paper (Sniegowska et al., in preparation) we will study a much
broader parameter range. However, they show the critical role of the outer
disk radius, and, on the other hand, a moderate dependence of the
modifications of the description of the disk structure. This second aspect is
very promising for the future, since the disk structure, corona formation and
the transition to ADAF are still subject of vigorous modeling, both numerical
and analytical or semi-analytical.

Summarizing, we see from the simulations presented in this paper that the
timescales of the order of years can be recovered from the model, after proper
adjustment of the parameters. However, rapid oscillations discovered recently
in two sources and named the Quasi-Periodic Ejections seem to be beyond the
reach of the radiation pressure instability model. Below we include a short
discussion of these observed phenomena.

###  4.1 Case study - timescales in NGC 1566

NGC 1566 is one of the sources in which outbursts seem to appear semi-
repetitively, in the timescales of years, it is thus a source which may be
modeled by radiation pressure instability, as postulated by Sniegowska .
(2020). First observations of the source which report optical variability are
obtained by Quintana . (1975). Later, changes in H$\beta$ and H$\alpha$ lines
were detected, as well as the change in spectral continuum Alloin . (1985).
Monitoring of those changes was continued by Alloin . (1986).

New changes (hard X-ray emission) in this source were detected by Ducci .
(2018) with INTEGRAL (Winkler ., 2003) mission and then followed up with using
SWIFT observatory data by Ferrigno . (2018). Those observations were
complemented by photometry and spectroscopy later on. Oknyansky . (2019)
report the increase in the luminosity by 25-30 times in June 2018, which was
followed by smaller outbursts and makes this source even more interesting.

Jana . (2021) analyse XMM-Newton (Jansen ., 2001), NuSTAR (Harrison ., 2013)
and Swift (Evans ., 2009) X-ray band archival observations of NGC 1566 from
the period between 2015 and 2019.

Amplitudes in the X-ray band change by the factor of 70, whereas bolometric
luminosity by the factor of 25. Our models shown in Fig. 4 imply higher
amplitude in the UV band, but searching the parameter space we can possibly
find more suitable solutions. The timescales of the order of a few years are
also longer than what we show in Fig. 4, but allowing for higher viscosity
parameter, adjusting the parameters in such a way as to get lower outburst
amplitude, and (more importantly) for still lower outer radius we could
shortened the period. Using the timescale for radiation pressure dominated
outburst from Grzȩdzielski . (2017) for the 107 solar masses and change in
luminosity factor as 25 for this source we get 135 years as the minimal
timescale of outburst. However, in this paper modified viscosity was used
(which in general shortened the period) but the outer radius was large
allowing for unconstrained radial development of the outburst. A smaller outer
radius in these models also could easily shorten the period to the requested
few years timescale.

###  4.2 Case study - timescales in RX J1301.9+2747

The regular quasi-periodic eruptions were discovered so far in two sources
(Miniutti ., 2019; Giustini ., 2020). The second of these two papers brings
the discovery of regular X-ray outbursts in the timescales of the order of 5
hours. The black hole mass in this source is lower than the value used in our
simulations ($\sim 10^{6}M_{\odot}$ vs. $\sim 10^{7}M_{\odot}$, so the
timescales from the model should be scaled down by (at least) the mass ratio
but such short timescales cannot be reached within the frame of the radiation
pressure instability model considered in this paper. As discussed by Giustini
. (2020), the underlying mechanism requires a large height to radius ratio,
i.e. the instability must happen not in the cold disk zone but (most likely)
in the warm corona. In our model, we have only hot corona at local virial
temperature which is stable, and the warm corona is not present. Hints for
thermal instabilities come from the problems to derive stationary warm corona
solutions for some parameter range (see Gronkiewicz Różańska, 2020) but no
time-dependent models are available so far.

##  5 Summary

In this work:

  * •

we perform the full computations of time-dependent evolution of the accretion
disk induced by the radiation pressure using GLADIS code,

  * •

we confirm that it is possible to obtain limit-cycle oscillations for the
black hole mass 107M⊙ using this model,

  * •

inclusion of the inner ADAF and coronal flow does not shorten the outburst,

  * •

the timescales from the model are not very sensitive to the disk/corona
coupling if outbursts are present, but the amplitudes do change,

  * •

the outburst timescales shorten strongly with the decrease of the outer disk
radius,

  * •

the phenomenon of CL AGN seems to be complex and challenging to model, since
observed objects show different behaviors. The radiation pressure model may
describe (after further modifications) the yearly variability in the sources
like NGC 1566, but most probably a different mechanism must account for Quasi-
Periodic Ejections.

However, these results require confirmation, with testing the different
resolution. The code is grid-based, and such hydrodynamical codes are known to
have resolution-dependent convergence.

## Acknowledgments

MS acknowledges the organizers and participants of 13th Serbian Conference on
Spectral Line Shapes in Astrophysics for the fruitful discussion and comments
which help to improve this work. BC and MS acknowledge the financial support
from the Polish Funding Agency National Science Centre, project
2017/26/A/ST9/-00756 (MAESTRO 9) and MNiSW grant DIR/WK/2018/12.

The package (tar.gz) with the source code GLADIS is available at
$http://www.cft.edu.pl/astrofizyka/$ under tab Numerics.

## References

  * Alloin . (1985) 1985ApJ…288..205AAlloin, D., Pelat, D., Phillips, M. Whittle, M.  198501,  ApJ288205-220.  10.1086/162783
  * Alloin . (1986) 1986ApJ…308…23AAlloin, D., Pelat, D., Phillips, MM., Fosbury, RAE. Freeman, K.  198609,  ApJ30823.  10.1086/164475
  * Czerny (2019) 2019Univ….5..131CCzerny, B.  201905,  Universe55131.  10.3390/universe5050131
  * Ducci . (2018) 2018ATel11754….1DDucci, L., Siegert, T., Diehl, R., Sanchez-Fernandez, C., Ferrigno, C., Savchenko, V. Bozzo, E.  201806,  The Astronomer’s Telegram117541.
  * Evans . (2009) 10.1111/j.1365-2966.2009.14913.xEvans, PA., Beardmore, AP., Page, KL. et al.  200907,  Monthly Notices of the Royal Astronomical Society39731177-1201.  https://doi.org/10.1111/j.1365-2966.2009.14913.x 10.1111/j.1365-2966.2009.14913.x
  * Feng . (2021) 2021ApJ…916…61FFeng, J., Cao, X., Li, Jw. Gu, WM.  202108,  ApJ916261.  10.3847/1538-4357/ac07a6
  * Ferrigno . (2018) 2018ATel11783….1FFerrigno, C., Siegert, T., Sanchez-Fernandez, C., Kuulkers, E., Ducci, L., Savchenko, V. Bozzo, E.  201806,  The Astronomer’s Telegram117831.
  * Giustini . (2020) Giustini2021Giustini, M., Miniutti, G. Saxton, RD.  202004,  A&A636L2.  10.1051/0004-6361/202037610
  * Graham . (2020) 2020MNRAS.491.4925GGraham, MJ., Ross, NP., Stern, D. et al.  202002,  MNRAS49144925-4948.  10.1093/mnras/stz3244
  * Gronkiewicz Różańska (2020) Gronkiewicz2020Gronkiewicz, D. Różańska, A.  202001,  A&A633A35.  10.1051/0004-6361/201935033
  * Grzȩdzielski . (2017) grzedzielski2017Grzȩdzielski, M., Janiuk, A., Czerny, B. Wu, Q.  201707,  A&A603A110.  10.1051/0004-6361/201629672
  * Harrison . (2013) 2013ApJ…770..103HHarrison, FA., Craig, WW., Christensen, FE. et al.  201306,  ApJ7702103.  10.1088/0004-637X/770/2/103
  * Hutsemékers . (2019) 2019A&A…625A..54HHutsemékers, D., Agís González, B., Marin, F., Sluse, D., Ramos Almeida, C. Acosta Pulido, JA.  201905,  A&A625A54.  10.1051/0004-6361/201834633
  * Jana . (2021) 2021MNRAS.507..687JJana, A., Kumari, N., Nandi, P. et al.  202110,  MNRAS5071687-703.  10.1093/mnras/stab2155
  * Janiuk (2020) 2020mbhe.confE..48JJaniuk, A.  202012,  GLADIS: GLobal Accretion Disk Instability Simulation GLADIS: GLobal Accretion Disk Instability Simulation.  Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. Palermo Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. Palermo 48.
  * Janiuk Czerny (2005) 2005MNRAS.356..205JJaniuk, A. Czerny, B.  200501,  MNRAS3561205-216.  10.1111/j.1365-2966.2004.08435.x
  * Janiuk Czerny (2007) 2007A&A…466..793JJaniuk, A. Czerny, B.  200705,  A&A4663793-803.  10.1051/0004-6361:20066133
  * Janiuk Czerny (2011) 2011MNRAS.414.2186JJaniuk, A. Czerny, B.  201107,  MNRAS41432186-2194.  10.1111/j.1365-2966.2011.18544.x
  * Janiuk . (2002) janiuk2002Janiuk, A., Czerny, B. Siemiginowska, A.  200209,  ApJ576908-922.  10.1086/341804
  * Janiuk . (2015) 2015janiukJaniuk, A., Grzedzielski, M., Capitanio, F. Bianchi, S.  201502,  A&A574A92.  10.1051/0004-6361/201425003
  * Jansen . (2001) 2001A&A…365L…1JJansen, F., Lumb, D., Altieri, B. et al.  200101,  A&A365L1-L6.  10.1051/0004-6361:20000036
  * Kokubo Minezaki (2020) 2020MNRAS.491.4615KKokubo, M. Minezaki, T.  202002,  MNRAS49144615-4633.  10.1093/mnras/stz3397
  * Kozłowski . (2016) kozlowski2016Kozłowski, S., Kochanek, CS., Ashby, MLN. et al.  2016Feb,  ApJ817119.  10.3847/0004-637X/817/2/119
  * Lawrence Papadakis (1993) lawrence1993Lawrence, A. Papadakis, I.  1993Sep,  ApJ414L85.  10.1086/187002
  * Lightman Eardley (1974) le74Lightman, AP. Eardley, DM.  197401,  ApJ187L1.  10.1086/181377
  * MacLeod . (2016) 2016MNRAS.457..389MMacLeod, CL., Ross, NP., Lawrence, A. et al.  201603,  MNRAS4571389-404.  10.1093/mnras/stv2997
  * Miniutti . (2019) Miniutti2019Miniutti, G., Saxton, RD., Giustini, M. et al.  201909,  Nature5737774381-384.  10.1038/s41586-019-1556-x
  * Noda Done (2018) noda2018Noda, H. Done, C.  2018Nov,  MNRAS4803898-3906.  10.1093/mnras/sty2032
  * Oknyansky . (2019) 2019MNRAS.483..558OOknyansky, VL., Winkler, H., Tsygankov, SS. et al.  201902,  MNRAS4831558-564.  10.1093/mnras/sty3133
  * Pringle . (1973) 1973A&A….29..179PPringle, JE., Rees, MJ. Pacholczyk, AG.  197312,  A&A29179.
  * Quintana . (1975) 1975MNRAS.173P..57QQuintana, H., Kaufmann, P. Sasic, JL.  197511,  MNRAS17357P-60P.  10.1093/mnras/173.1.57P
  * Sánchez-Sáez . (2021) 2021arXiv210607660SSánchez-Sáez, P., Lira, H., Martí, L. et al.  202106,  arXiv e-printsarXiv:2106.07660.
  * Scepi . (2021) 2021MNRAS.502L..50SScepi, N., Begelman, MC. Dexter, J.  202103,  MNRAS5021L50-L54.  10.1093/mnrasl/slab002
  * Sesar . (2007) sesar2007Sesar, B., Ivezić, Ž., Lupton, RH. et al.  200712,  AJ1342236-2251.  10.1086/521819
  * Shakura Sunyaev (1973) ss73Shakura, NI. Sunyaev, RA.  1973,  A&A24337-355.
  * Sniegowska . (2020) 2020A&A…641A.167SSniegowska, M., Czerny, B., Bon, E. Bon, N.  202009,  A&A641A167.
  * Wang Bon (2020) 2020wangbonWang, JM. Bon, E.  202011,  A&A643L9.  10.1051/0004-6361/202039368
  * Winkler . (2003) 2003A&A…411L…1WWinkler, C., Courvoisier, TJL., Di Cocco, G. et al.  200311,  A&A411L1-L6.  10.1051/0004-6361:20031288
  * Woo Urry (2002) woo2002Woo, JH. Urry, CM.  200211,  ApJ579530-544.  10.1086/342878
  * Zhang (2021) 2021MNRAS.500L..57ZZhang, XG.  202101,  MNRAS5001L57-L61.  10.1093/mnrasl/slaa184


# Kilo-Tesla axial magnetic field generation with high intensity spin and
orbital angular momentum beams

Andrew Longman  Lawrence Livermore National Laboratory, Livermore, California
94551<EMAIL_ADDRESS>Robert Fedosejevs  Department of Electrical and
Computer Engineering, University of Alberta, Edmonton T6G1R1, Canada

###### Abstract

Absorption of angular momentum from a high intensity laser pulse can lead to
the generation of strong axial magnetic fields in plasma. The effect, known as
the inverse Faraday effect can generate kilo-Tesla strength, multi-picosecond,
axial magnetic fields extending over hundreds of microns in underdense plasma.
In this paper we explore the effect with ultra-high intensity circularly
polarized Gaussian beams and linearly polarized orbital angular momentum beams
comparing analytic expressions with 3D particle-in-cell simulations. We
develop a model for the transverse magnetic field profiles, introduce a new
model for the temporal decay, and show that while the magnetic field strength
is independent of plasma density, it has a strong dependence on the laser beam
waist.

##  I Introduction

The inverse Faraday effect (IFE) describes the generation of axial magnetic
fields when angular momentum is transferred from a laser pulse to plasma. If
the pulse is circularly polarized (CP) and of relativistic intensity, the
generated magnetic field strength can be on the order of 100’s of Tesla, last
several picoseconds, and extend over millimeter scales [1]. The IFE has been
studied extensively with CP beams both theoretically [2, 3, 4, 5, 6, 7, 8],
and experimentally [1, 9, 10]; yet the effect is not exclusive to CP beams but
also applicable to lasers carrying orbital angular momentum (OAM) [11, 12,
13]. Currently, we only know of one analytic model describing OAM driven
magnetic fields and its transverse profile [11], but has not been verified
numerically or experimentally. Furthermore, little is known about the axial
extent of the magnetic fields or their lifetime, for both CP and OAM drivers.

There are several advantages of using OAM beams in high powered lasers as
opposed to CP beams: unlike CP beams that carry a maximum spin angular
momentum of $\sigma_{z}\hbar$ per photon where $\sigma_{z}=\\{\pm 1\\}$ is the
spin number ($\sigma_{z}=0$ for linearly polarized (LP) beams), OAM beams
carry an angular momentum of $(\sigma_{z}+\ell)\hbar$ per photon where $\ell$
is an azimuthal mode number that can take any integer value [14]. OAM has also
been shown to couple more efficiently to free electrons due to an optimal
overlap between the laser intensity profile and the angular momentum density
yielding greater coupling at larger radial distances [15]. Additionally, high
intensity OAM beams have recently been produced in high power lasers using
off-axis spiral phase mirrors mitigating nonlinear effects like temporal
distortions and B-integral present with transmissive optics such as quarter
waveplates [16].

Recent works have looked at simulating IFE driven magnetic fields from OAM
beams in various configurations: OAM beams with radial and azimuthal
polarization [12], amplification of seeded magnetic fields [17],
spatiotemporal light springs [18], and most recently linearly polarized OAM
beams [13]. In these studies the laser intensities were of moderate intensity
$(I_{0}\approx 10^{18}$ Wcm${}^{-2})$ verifying the existence of weaker
magnetic fields $(\approx 10$ T$)$, with little modelling of the spatial or
temporal properties of the magnetic fields.

In this work, we explore the spatial, and for the first time the temporal
scaling of IFE magnetic fields driven by ultra-high intensity $(I_{0}\approx
10^{20}$ Wcm${}^{-2})$ CP Gaussian and LP OAM beams. We develop new analytic
models for both the spatial and temporal scales of the magnetic field and
verify them with full 3D relativistic particle-in-cell (PIC) simulations. We
demonstrate strong coupling of laser angular momentum to plasma through
ponderomotive forces, resulting in axial magnetic fields in excess of 1 kT,
more than 200 $\mu$m in length, and lasting several picoseconds in good
agreement with the predictions of our analytic models, and previous
experiments using CP Gaussian beams.

##  II Coupling of high intensity angular momentum to plasma

A circularly polarized beam is classically defined by its polarization vector
rotating about the laser propagation axis $k_{z}$ in either a clockwise or
counter-clockwise direction. Similarly, a beam carrying a well defined OAM is
typically characterized by a helical wavefront also orientated about the
$k_{z}$ vector with the helix rotating in either a clockwise or counter-
clockwise direction [14]. Unlike a Gaussian beam, paraxial focusing of an OAM
beam produces a complicated field structure described by modified-Bessel or
hypergeometric functions but can be modelled using a suitably chosen basis set
such as the Laguerre-Gaussian (LG) modes [19]. The electric field magnitude of
a general OAM laser pulse can then be approximately given by
$E=E_{0}\psi_{\ell}g(t)\exp[i(\omega_{0}t-k_{z}z)]$ where $E_{0}$ is the peak
electric field strength of the fundamental Gaussian mode ($\ell=0$), $g(t)$ is
the temporal envelope, and $\omega_{0}$ is the laser frequency. $\psi_{\ell}$
are the azimuthal LG modes (zero-radial mode) and can be given in the
following compact form [14, 20],

$\psi_{\ell}=\frac{1}{\sqrt{|\ell|!}}\left(\frac{r\sqrt{2}}{w_{0}}\right)^{|\ell|}\left(\frac{z_{0}}{Z}\right)^{|\ell|+1}exp\left(-\frac{kr^{2}}{2Z}\right)e^{i\ell\phi},$
(1)

where $w_{0}$ is the beam waist in the focal plane, $z_{0}=kw_{0}^{2}/2$ is
the Rayleigh range, and $Z=z_{0}+iz$ is the complex beam parameter. Taking the
modulus squared of the field amplitude given in Eq. 1 for $|\ell|\geq 1$, we
find the so called “donut” mode intensity profiles that are symmetric about
the $k_{z}$ axis; the helical phase is lost in the modulus resulting in no
azimuthal structure in the intensity profile [14]. Increasing the azimuthal
mode integer $|\ell|$ increases the OAM in the beam, increases the donut
radius and by conservation of energy, the peak intensity around the donut
decreases for a fixed laser energy and focusing geometry [21]. As the peak
intensity varies for different values of $|\ell|$ in a non-trivial way, we opt
to use the beam power instead as it is constant for all $|\ell|$ modes.

The ponderomotive force acting on an electron associated with a high intensity
focal spot can be given by,
$\bm{F}_{p}=-m_{e}c^{2}\bm{\nabla}\sqrt{1+|\bm{a}|^{2}}$ where
$\bm{a}=e\bm{E}/m_{e}c\omega_{0}$ is the normalized vector potential, $e$ and
$m_{e}$ are the electron charge and mass respectively, and $c$ is the speed of
light. In the standard case of a high-intensity Gaussian beam a finite energy
maybe transferred from the laser to plasma electrons through ponderomotive
scattering [22]. The same is true for higher order LG modes [23], with some
possibility for electron trapping and additional energy transfer within the
donut mode leading to a possible increased absorption rate [24, 25]. However,
ponderomotive scattering to first order is an intensity driven phenomenon, and
no OAM is transferred to the plasma as the LG mode intensity has no azimuthal
component.

Recent work has analyzed the motion of free electrons in a LG beam in more
detail, expanding to first and second order perturbations [13, 15]. It was
shown that angular momentum can be transferred to free electrons in LP LG
beams when considering second order terms in the equations of motion; the
electron gains and loses angular momentum but averaging over one laser cycle
yields zero net angular momentum transfer. An asymmetry in the laser field
such as the ponderomotive force, or collisions are therefore needed to
transfer a net OAM to the electrons.

##  III The Inverse Faraday Effect

To model the IFE, we start with the angular momentum conservation model used
in Haines [3]. In this model, coupling of the angular momentum density in the
plasma with that of a driving laser is performed through an arbitrary
absorption parameter yielding the conservation equation,

$\displaystyle n_{e}m_{e}r\frac{dv_{e\theta}}{dt}=-n_{e}erE_{\theta}$
$\displaystyle-
n_{e}er(v_{ez}B_{r}-v_{er}B_{z})+\alpha_{abs}M_{z}c-n_{e}m_{e}\nu_{ei}rv_{e\theta}.$
(2)

Here, $n_{e}$ is the electron density, $\alpha_{abs}$ is the laser absorption
fraction per unit length, $M_{z}$ is the laser angular momentum density,
$v_{er},v_{e\theta}$ are the radial and azimuthal electron velocities
respectively, and $\nu_{ei}$ is the electron-ion collision frequency.
$E_{\theta},B_{r},$ and $B_{z}$ are the plasma azimuthal electric, radial
magnetic and axial magnetic fields respectively. Assuming the electrons are in
steady state, ignoring the inertial and collisional terms, and using Faraday’s
law, we obtain the time rate of change of the axial magnetic field
approximately as [3],

$\frac{\partial B_{z}}{\partial t}\approx\frac{c}{re}\frac{\partial}{\partial
r}\frac{\alpha_{abs}M_{z}}{n_{e}},$ (3)

For an OAM mode that is either CP or LP, the angular momentum density can be
given by [11, 14],

$M_{z}=\frac{I_{0}|\psi_{\ell}|^{2}g(t)^{2}}{\omega_{0}c}\left[\ell+\sigma_{z}\left(|\ell|-\frac{2r^{2}}{w_{0}^{2}}\right)\right].$
(4)

Radial differentiation of Eq. 4 is straightforward, and if we assume a
Gaussian temporal function, then the integration of Eq. 3 is trivial given
$\int g(t)^{2}dt\approx 3\tau/4$ where $\tau$ is the temporal full-width half
maximum (FWHM). By introducing the pulse temporal function $g(t)$, we are
limiting the IFE to only be driven when the laser field is present. Though the
magnetic field does not instantaneously disappear after the laser has passed,
its evolution is no longer governed by the IFE, but rather the plasma dynamics
which will be discussed later.

Substituting Eq. 1 into Eq’s. 3 and 4, assuming the laser absorption rate and
electron density is constant in space and time (no ponderomotive channelling),
a Gaussian temporal pulse shape, and substituting the laser power $P=\pi
I_{0}w_{0}^{2}/2$ we obtain the following expression for the axial magnetic
field,

$\displaystyle B_{z}=\frac{3\alpha_{abs}P\tau}{\pi
en_{e}\omega_{0}w_{0}^{4}}|\psi_{\ell}|^{2}\times$
$\displaystyle\left[\ell\left(\frac{|\ell|w_{0}^{2}}{r^{2}}-2\right)+\sigma_{z}\left(\frac{|\ell|^{2}w_{0}^{2}}{r^{2}}-4|\ell|+\frac{4r^{2}}{w_{0}^{2}}-2\right)\right].$
(5)

The direct proportionality on $\ell$ and $\sigma_{z}$ in Eq. III indicates
that the axial magnetic field direction can be controlled by changing the sign
of OAM helicity, or the CP handedness. The helicity of a LP OAM mode, the
magnetic field direction and the corresponding electron motion is illustrated
in Fig. 1.

Figure 1:  Illustration of an IFE driven B-field with a LP $\ell=1$ OAM mode
and the corresponding electron helical trajectory in red with an arrow head
showing the velocity direction. The green and yellow isosurfaces show the
positive and negative electric fields of the laser respectively while the blue
isosurface shows the positive magnetic field direction as indicated by the
arrow.

Several transverse magnetic field profiles are plotted in Fig. 2 with laser
plasma parameters: $n_{e}=0.03n_{c}$, $\lambda=1$ $\mu$m, $w_{0}=6$ $\mu$m,
$\alpha_{abs}=0.435$ mm-1, $\tau$=100 fs, and $P=$ 65 TW, where $\lambda$ is
the laser wavelength, and $n_{c}$ is the critical plasma density. We note that
the $(|\sigma_{z}|=1,|\ell|=0)$ and $(\sigma_{z}=0,|\ell|=1)$ modes produce
the same magnetic field profiles peaked on axis, while the higher order $\ell$
modes produce more complex coaxial structures that are zero on axis. Given the
magnetic field peak on axis and the higher magnetic field strength, the
$(|\sigma_{z}|=1,|\ell|=0)$ and $(\sigma_{z}=0,|\ell|=1)$ driven fields are
probably more desirable for most applications. For the purpose of this work,
we do not consider beams that are both circularly polarized and have OAM, but
we remark that a beam with $(\sigma_{z}=\pm 1,\ell=\pm 1)$ produces the
strongest on-axis magnetic field with a peak field strength double that of the
$(|\sigma_{z}|=1,|\ell|=0)$ case. Using these parameters, and assuming the
laser is entirely absorbed, the peak ratio of the energy per unit length
contained in the magnetic field to the laser energy per unit length is roughly
$0.2\%$ for the CP $\ell=0$ and LP $|\ell|=1$ cases, slightly increasing for
higher order $\ell$ modes.

Figure 2:  Transverse IFE driven magnetic field profiles driven from a CP
Gaussian beam and various LP OAM beams given by Eq. III and the legend. The
following laser plasma parameters were used: $n_{e}$=0.03$n_{c}$, $\lambda=1$
$\mu$m, $P=65$ TW, $w_{0}=6$ $\mu$m, $\alpha_{abs}=0.435$ mm-1, and $\tau=100$
fs.

The dependence of $\alpha_{abs}$ on $\ell$ is not currently understood for any
absorption mechanism. At laser intensities of $I\lambda^{2}\geq 10^{19}$
Wcm${}^{-2}\mu$m2 and pulse duration’s greater than 100 fs, we can estimate
the pump depletion length in underdense plasma ($n_{e}>0.01n_{c}$) as
$L_{pd}=n_{c}c\tau/n_{e}$ where $\tau$ is the temporal FWHM [26, 27].
Rearranging, we obtain the absorption rate as,

$\alpha_{abs}=\frac{n_{e}}{2n_{c}c\tau}.$ (6)

Given the absorption is independent of beam intensity and beam waist, we
assume the model is approximately valid for low order OAM modes. If
$\tau\gtrsim 100$ fs, and $n_{e}\gtrsim 0.01n_{c}$, we can use Eq’s. III and 6
to estimate the peak axial magnetic field strength of the
$(|\sigma_{z}|=1,|\ell|=0)$ and $(\sigma_{z}=0,|\ell|=1)$ modes,

$|B|_{max}\approx 10\frac{P[TW]\lambda^{3}[\mu m]}{w_{0}^{4}[\mu m]}kT.$ (7)

The fourth power dependence on the laser beam waist indicates a strong
dependence on the beam f-number and relativistic self-focusing. Eq. 7 should
therefore be considered an under-estimate as the beam will self-focus beyond
the diffraction limit giving rise to much stronger magnetic fields. Using Eq.
7, we can compare its prediction with the experimental results previously
measured using CP Gaussian beams [1]. In that experiment, the Vulcan laser
($40$ TW, $1$ ps, $\lambda=1.054\mu$m) was focussed into an underdense plasma
where it was observed to self focus to a beam waist of $w_{0}\approx 5\mu$m. A
peak magnetic field strength was observed of $\approx 700$ T, in good
agreement with our Eq. 7 estimate of $771$ T. The same experiment also
measured an independence of the magnetic field strength on the plasma density,
also in agreement with our model.

Circularly polarizing large diameter, high power laser beams is not trivial
requiring fragile and costly optics. However, with the new generation of PW
class lasers coupled with off-axis spiral phase mirrors [16], axial magnetic
fields on the order of 10’s of kT may be feasible using LP OAM modes.
Ultimately, the field strength will depend on the absorption rate of the laser
which may diverge from our simple model when using ultra-high intensity beams
$(I_{0}\gtrsim 10^{21}Wcm^{-2})$, high $|\ell|$ beams, and ultra-short pulse
durations $(\tau<100fs)$. Additionally, realistic diffraction of OAM beams
when generated by spiral phase optics is more complex leading to larger beam
waists and lower intensities than one would expect using pure Gaussian beams
[21, 19]. While we do not include these modified beam waists and intensities
here, it is simple to use the estimates from these previous works to estimate
the reduction in the magnetic field strength.

Figure 3:  Simulation results for CP Gaussian driven magnetic fields at times
(from left to right) $0.6$ ps, $0.77$ ps, and $0.93$ ps. Tiles a), d), and g):
longitudinal electron density slice maps normalized to the initial electron
density. Tiles b), e), and h): longitudinal axial magnetic field slice maps
averaged over 33 fs, the laser pulse is overlayed in green and lineout regions
are given by the red dashed lines. Tiles c), f), and i): Radial axial magnetic
field lineouts averaged azimuthally, temporally, and longitudinally between
the red dashed lines in the corresponding tile above, given by the red solid
lines. The theoretical prediction of the model given by Eq. III are shown by
blue dashed lines. Laser-plasma parameters are given in the main text.

##  IV Simulation results

To simulate the IFE, we use the 3D relativistic PIC code EPOCH [28]. We assume
a fully ionized helium plasma with a super-Gaussian longitudinal shape to
mimic the electron density from a gas jet given by,
$n(z)=n_{e}\exp\\{-[(z-350\mu m)/300\mu m]^{10}\\}$ with initial electron and
ion temperatures of 1 keV and 1 eV respectively. Using this profile, the
plasma ramps from vacuum to $n_{e}$ over roughly 100 $\mu$m. The laser is
polarized along the $\hat{x}$ axis, has a Gaussian temporal pulse shape, and
is focused to the plane $z=135$ $\mu$m.

The laser pulse is tracked with a moving window 200 $\mu$m long and 48 $\mu$m
$\times$ 48 $\mu$m in the transverse directions with grid cell sizes of 50 nm
$\times$ 80 nm $\times$ 80 nm respectively. We use 4 particles per cell and
have open boundaries throughout. After the magnetic field has been generated
and before significant laser diffraction and filamentation occurs, the moving
window stops at 1 ps and the evolution of the magnetic field is observed.
Simulations were run for the $(\sigma_{z}=0,\ell=0)$, $(\sigma_{z}=1,\ell=0)$,
$(\sigma_{z}=0,\ell=1)$ and $(\sigma_{z}=0,\ell=2)$ laser modes, with
parameters $P=65$ TW, $\tau=100$ fs, $\lambda=1$ $\mu$m, $w_{0}=6$ $\mu$m, and
$n_{e}=3\times 10^{19}$ cm-3. A fifth simulation with $(\sigma_{z}=0,\ell=1)$,
$P=50$ TW, $\tau=100$ fs, $\lambda=1$ $\mu$m, $w_{0}=7.2$ $\mu$m, and
$n_{e}=3\times 10^{19}$ cm-3 was also run to verify the decay model.

Figure 4:  Simulation results for LP $\ell=1$ Laguerre Gaussian driven
magnetic fields at times (from left to right) $0.6$ ps, $0.77$ ps, and $0.93$
ps. Tiles a), d), and g): longitudinal electron density slice maps normalized
to the initial electron density. Tiles b), e), and h): longitudinal axial
magnetic field slice maps averaged over 33 fs, the laser pulse is overlayed in
green and lineout regions are given by the red dashed lines. Tiles c), f), and
i): Radial axial magnetic field lineouts averaged azimuthally, temporally, and
longitudinally between the red dashed lines in the corresponding tile above,
given by the red solid lines. The theoretical prediction of the model given by
Eq. III are shown by blue dashed lines. Laser-plasma parameters are given in
the main text.  Figure 5:  Simulation results for LP $\ell=2$ Laguerre
Gaussian driven magnetic fields at times (from left to right) $0.6$ ps, $0.77$
ps, and $0.93$ ps. Tiles a), d), and g): longitudinal electron density slice
maps normalized to the initial electron density. Tiles b), e), and h):
longitudinal axial magnetic field slice maps averaged over 33 fs, the laser
pulse is overlayed in green and lineout regions are given by the red dashed
lines. Tiles c), f), and i): Radial axial magnetic field lineouts averaged
azimuthally, temporally, and longitudinally between the red dashed lines in
the corresponding tile above, given by the red solid lines. The theoretical
prediction of the model given by Eq. III are shown by blue dashed lines.
Laser-plasma parameters are given in the main text.

###  IV.1 Circularly Polarized Gaussian Simulations

A simulation was first run for the linearly polarized Gaussian beam to verify
the null magnetic field result, and also to verify the absorption model. While
we do not show the simulation results in this work the null field was
verified, and the absorption rate of this beam was found to be 0.43 mm-1 in
excellent agreement with the predicted value of 0.435 mm-1 from Eq. 6. A
second identical simulation was run but with a circularly polarized Gaussian
beam ($\sigma_{z}=1,\ell=0$) that measured absorption slightly enhanced to
0.48 mm-1.

Given the peak intensity of the CP Gaussian beam is $I_{0}=1.15\times
10^{20}Wcm^{-2}$, the interaction with the plasma is strongly nonlinear with
effects such as relativistic self-focussing, ponderomotive channeling, and
filamentation strongly dominating the interaction. This can make it
challenging to get clean magnetic field profiles out of the simulation,
especially at later times when the boundary conditions start to have an effect
on the simulation. Because of this, we examine the magnetic field at early
times before these nonlinear effects have time to sufficiently grow.

Fig. 3 shows the axial magnetic field generation from the CP Gaussian mode at
three early times in the simulation. Tiles a), d) and g) give the electron
density normalized to the initial electron density at times $0.6$ ps, $0.77$
ps, and $0.93$ ps respectively. Below these tiles are transverse slices of the
time-averaged axial magnetic fields (averaged over 33fs output dumps) relative
to the laser pulse overlayed in green for each of the three times. The bottom
three tiles give radial lineouts of the time averaged magnetic field that have
been averaged azimuthally, temporally (over 33fs), and longitudinally between
the red dashed lines in tiles b), e), and h). The simulation lineouts are
given in red, whereas the theoretical prediction of Eq. III is given with the
blue dashed line.

We find a remarkable overlap between the theory and the simulations at early
times but later diverges from the theory as the laser begins to diffract,
self-focus, and filament. In particular we see the peak magnetic field
strength enhance beyond the vacuum theory at $t=0.93$ ps due to the self-
focussing. We also note the persistence of the magnetic fields after the laser
pulse has passed due to the residual magnetization of the plasma.

###  IV.2 Linearly Polarized OAM Simulations

Next, simulations with LP OAM modes were run using both $\ell=1$ and $\ell=2$
modes, the results of which are given in Fig. 4 and Fig. 5 respectively with
the same layout described for the CP Gaussian beams in Fig. 3. The absorption
of the $\ell=1$ and $\ell=2$ modes were found to be 0.59 mm-1 and 0.60 mm-1
respectively. This increase in absorption may come from the larger cross
section of the OAM modes, electron trapping within the donut mode leading to a
more efficient laser-plasma energy transfer, or from a reduced ponderomotive
channeling leading to a higher plasma density interacting and coupling with
the laser directly.

Examining Fig. 4 we see a good overlap of the simulated magnetic field and
that predicted by the theory. Similar to the CP Gaussian simulations, we find
the early times of the $\ell=1$ simulation agree well, but at later times as
the laser starts to diffract, self-focus, and filament. In tile 4g) and h), we
see the magnetic field expand radially outwards without compromising its peak
strength. This surprising result could be related to the compressed electron
density on axis giving rise to a higher azimuthal current, and a slow radial
expansion outwards to fill the voids left by the ponderomotive force.

For the $\ell=2$ simulation shown in Fig. 5 we also see a good overlap between
theory and simulation at early times, but at later times we find that the
magnetic field amplifies to twice that of the predicted field, whilst still
maintaining its shape. This is likely due to the more complex self-focussing
of the LG modes as the donut mode both pinches into a tighter ring increasing
intensity but maintaining the beam waist, and then collapsing to a smaller
beam wasit at later times. During this self-focussing, additional angular
momentum may be coupled to the plasma giving rise to the enhanced magnetic
field. As our model does not include self-focussing, we do not expect to
reproduce the results at later times. There is some infilling of the null
within the magnetic field on axis as the plasma thermalizes, and then eventual
radial expansion as the laser later diffracts.

In the current configuration the plasma is radially isotropic and the axial
magnetic field length is essentially limited to the Rayleigh length of the
laser. For longer magnetic fields, one could use a preformed plasma channel to
guide the laser and mitigate the effects of diffraction but is not explored in
the current work.

Figure 6:  2D longitudinal slices of the time averaged (averaged over 33fs
output dumps) axial magnetic field and the corresponding time averaged
azimuthal electron velocity for the $(\sigma_{z}=0,\ell=1)$ OAM mode for times
$t=1,1.5,2$ ps. The green dashed lines indicate the regions used for sampling
the field decay in Fig. 7.

Fig. 6(a), (b) and (c) give the time averaged (averaged over 33fs output
dumps) axial magnetic field profile at times 1 ps, 1.5 ps and 2 ps
respectively for the $(\sigma_{z}=0,\ell=1)$ mode. We find the magnetic field
extends the length of the simulation box (200 $\mu$m) approximately equal to
the twice Rayleigh range ($z_{0}\approx 113$ $\mu$m).

At later times in Fig. 6(b) and (c), we find the $\ell=1$ mode to begin to
pinch, kink, and twist into a 3D spring-like shape. Similar kinking is found
for the $(|\sigma_{z}=1|,|\ell|=0)$ and $(|\sigma_{z}=0|,|\ell|=2)$ modes
although not shown. Fig 6 (d), (e), and (f) show transverse slices of the
electron azimuthal velocity time-averaged over 33 fs. We see the azimuthal
velocity map transition from a smooth cylindrical profile at 1 ps, to a more
turbulent and kinked map at 2 ps. This could be related to the wobble
instability found in theta-pinch configurations [29].

##  V Magnetic Field Decay

The lifetimes of the magnetic fields driven by the IFE at relativistic
intensities in collisionless plasmas are poorly understood. If we consider the
plasma beta, $(\beta=2\mu_{0}n_{e}k_{B}T/B^{2})$ where $\mu_{0}$ is the
permeability of free space, $k_{B}$ is Boltzmann’s constant, and $T$ is the
plasma temperature, we find that for magnetic field strengths of 1000 T,
electrons with energies $T<100$ keV are trapped in the magnetic field, while
higher energy electrons are able to escape leading to the subsequent decay of
the magnetic field.

To model the decay, we consider the motion of the hot electrons in helical
motion with both azimuthal and radial velocity components. From Fig. 6 d) -
f), we find the bulk of the hot electrons with OAM are born within a cylinder
of radius approximately equal to the peak intensity radius of the OAM mode
$r_{cyl}\approx w_{0}\sqrt{\ell/2}$. The higher order $|\ell|\geq 2$ modes are
more complex with coaxial fields and will not be considered in our model.

The magnetic field as a function of time from a rotating cylinder of plasma
can be derived from the Biot-Savart law [30]. The current density can be given
as $\bm{j}=qn_{e}(v_{r}\hat{r}+v_{\theta}\hat{\theta})$, the position vector
as $\bm{r}(t)=(v_{r}t+r_{0})\hat{r}+z\hat{k}$ where $r_{0}$ is the initial
electron radius, and $\hat{r}$ is the rotating polar coordinate. Inserting
these into the Biot-Savart equation we find,

$B(t)=\frac{\mu_{0}en_{e}}{2}\int_{0}^{c\tau}dz\int_{0}^{r_{cyl}}\frac{v_{\theta}(v_{r}t+r_{0})r_{0}dr_{0}}{\left[(v_{r}t+r_{0})^{2}+z^{2}\right]^{3/2}}.$
(8)

We note that the axial length of the plasma for integration is restricted to
approximately one pulse length $c\tau$. As the plasma is collisionless, we can
assume the angular momentum acquired by the electrons is in steady state, that
is $L_{z}=m_{e}v_{\theta}r_{0}$ is constant. Solving Eq. 8 is straightforward
if we perform the radial integral first yielding the result,

$B(t)=B_{0}\left[\sinh^{-1}\left(\frac{c\tau}{v_{r}t}\right)-\sinh^{-1}\left(\frac{c\tau}{v_{r}t+r_{cyl}}\right)\right].$
(9)

Here, $B_{0}$ is the peak axial magnetic field at genesis.

To estimate the radial velocity of the electrons we can use the gradient of
the plasma pressure [3]. To first order, the pressure can be attributed solely
to the thermal pressure given by,

$\frac{n_{e}m_{i}}{Z}\frac{dv_{r}}{dt}=-\frac{\partial}{\partial
r}\left(\frac{2}{3}\alpha_{abs}\int I(t)dt\right).$ (10)

Here, $m_{i}$ and $Z$ are the ion mass and ionization state respectively. By
using the ion mass, we assume the plasma to be quasineutral which was
previously shown to be approximately valid on these timescales and laser
intensities [3]. Integrating the laser temporal function is trivial as shown
earlier in the paper leaving a second order ODE,

$\frac{d^{2}r}{dt^{2}}+\frac{Z\alpha_{abs}I_{0}\tau}{2n_{e}m_{i}}\frac{\partial}{\partial
r}|\psi_{\ell}|^{2}=0.$ (11)

There are no known analytic solutions to Eq. 11 and we opt instead to use a
parabolic approximation of $\psi_{\ell}$ given by,

$|\psi_{\ell}|^{2}\approx
A-\left(\frac{r}{w_{0}}-\sqrt{\frac{\ell}{2}}\right)^{2},$ (12)

where $A$ is an arbitrary amplitude scaling variable that is lost when we use
the $\partial/\partial r$ operator.

The soltuion of the ODE is elementary using the parabolic approximation.
Averaging the radial velocity over all positions within a beam waist, we find

$\langle
v_{r}\rangle=\sqrt{\frac{Z\alpha_{abs}I_{0}\tau}{4\exp(2)n_{e}m_{i}}}.$ (13)

Substituting our previous simulation values, we find an average radial
velocity of $\langle v_{r}\rangle\approx 0.01c$. Given this, we can use the
approximation that $(c\tau)^{2}\gg(v_{r}t)^{2}$ given the timescale of the
magnetic field decay is observed to be on the order of a picosecond, allowing
us to simplify Eq. 9 to

$\frac{B(t-t_{0})}{B_{0}}\approx\ln\left(1+\frac{1}{t/t_{D}+0.582}\right).$
(14)

Here $t_{D}=w_{0}/(v_{er}\sqrt{2})$ is the decay parameter for the $|\ell|=1$
mode, and $t_{0}$ is the genesis time. The factor of 0.582 is set such that
$B(0)=B_{0}$. Combining the results of Eq. 13 and Eq. 14 allows us to derive
the scaling law for the time decay parameter of the $|\ell|=1$ driven magnetic
fields,

$t_{D}[ps]\approx\frac{1}{8}\frac{w_{0}^{2}[\mu m]}{\lambda[\mu
m]}\sqrt{\frac{A}{Z}}\frac{1}{\sqrt{P[TW]}},$ (15)

where $A$ is the atomic mass number of the plasma. While this model is not
valid for $|\ell|\geq 2$, it could be suitable for also modelling the CP
$\ell=0$ mode given its similarities to the $|\ell|=1$ magnetic field profile.

Using our simulation values, we estimate the decay parameter as $t_{D}=0.8$
ps. Determining the decay time numerically is challenging due to the magnetic
field instabilities, as well as laser self-focussing causing the magnetic
field to be dynamic on axis. We therefore radially, azimuthally, and
longitudinally average the magnetic field through a cylinder of radius
$r_{cyl}=w_{0}/\sqrt{2}$ and length $c\tau$. The longitudinal region is
indicated by the dashed green bars from $z=200-230$ $\mu$m as shown in Fig. 6.
The average value of the magnetic field in this region is plotted as a
function of time in Fig. 7 for three simulations. The simulation for the $65$
TW, $(\sigma_{z}=0,\ell=1),w_{0}=6$ $\mu$m mode is shown with blue circles
where as the $50$ TW, $(\sigma_{z}=0,\ell=1)$ and $w_{0}=7.2$ $\mu$m
simulation is shown with red diamonds. The purple crosses represent the data
from the initial CP Gaussian $(\sigma_{z}=1,\ell=0)$ simulation with $P=65$
TW, $\ell=0,w_{0}=6$ $\mu$m.

Figure 7:  Average magnetic field within a cylinder of radius
$r_{cyl}=w_{0}/\sqrt{2}$ and length 30 $\mu$m from $z-200-230\mu$m as a
function of time. The LP $65$ TW, $\ell=1,w_{0}=6$ $\mu$m mode is given by the
blue circles, and fitted with the red dashed line with decay parameter
$t_{D}=0.86$ ps. The LP $50$ TW, $\ell=1,w_{0}=7.2$ $\mu$m mode is given by
the red diamonds, and fitted with the green dashed line with decay parameter
$t_{D}=1.18$ ps. The $65$ TW, $\sigma_{z}=1,\ell=0,w_{0}=6$ $\mu$m mode is
given by the purple crosses, and fitted with the black dashed line with decay
parameter $t_{D}=0.69$ ps.

Using a nonlinear least squares fit of Eq. 14 to the data, we find a numerical
decay parameter of $0.86$ ps for the $\sigma_{z}=0,\ell=1,w_{0}=6$ $\mu$m
mode, shown by the red dashed line in Fig. 7. This is in very good agreement
with the analytic estimate of $0.8$ ps. Considering the
$\sigma_{z}=0,\ell=1,w_{0}=7.2$ $\mu$m case with a laser power of $50$ TW, we
predict a decay parameter of $\tau=1.3$ ps and numerically find a value of
$1.18$ ps shown by the green dashed line in Fig. 7.

Finally for the $(\sigma_{z}=1,\ell=0)$ case we find a predicted decay also of
$0.8$ ps, and numerically find a decay time of 0.69 ps as shown by the black
dashed line. This reduced decay time could be due to the stronger self-
focussing of the CP Gaussian mode as it has a higher peak intensity than the
LG modes. The increased self-focussing results in a smaller beam waist, and
hence a shorter decay time.

Given the dependence of the decay time on the laser beam waist, one could
control the magnetic field lifetime by adjusting the f-number to higher values
at a cost of magnetic field strength. An optimal f-number is therefore needed
to balance the magnetic field strength, axial length, and decay time to suit a
given application of interest. Alternatively one could adjust the decay time
of the magnetic field by changing the plasma to a heavier ion species.

##  VI Summary

In summary, we have successfully demonstrated and characterized the IFE with
CP Gaussian and linearly polarized OAM modes over large spatial and temporal
scales. We have derived a description of OAM driven magnetic fields in
underdense plasma and introduced a model for the subsequent decay of the
magnetic field. Magnetic fields with strengths up to 1 kT extending upto $200$
$\mu$m and persisting for several picoseconds have been demonstrated.
Simulations indicate an increased absorption rate when using increased OAM
mode numbers. The decay time, axial length, and magnetic field strength are
all functions of the laser beam waist and can be optimized to suit an
experiment as such. Plasma channel guiding of the laser pulse could extend the
magnetic field over several Rayleigh lengths and would instead be limited to
the pump depletion length.

With the recent demonstration of high intensity OAM modes in high-power laser
facilities with parameters similar to those assumed in this work, experimental
verification of the magnetic fields could now be feasible [16]. Generation and
control of multi-kT magnetic fields will open up new opportunities in a number
of areas of high energy density physics and laboratory astrophysics including
particle acceleration schemes, magnetic re-connection, fast ignition and
fundamental physics.

###### Acknowledgements.

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344,
supported by the Natural Sciences and Engineering Research Council of Canada
(RGPIN-2019-05013), and support provided by WestGrid and Compute Canada. This
work was in part funded by the UK EPSRC grants EP/G054950/1, EP/G056803/1,
EP/G055165/1 and EP/ M022463/1. The authors would like to thank and
acknowledge J. Myatt, J. Ludwig, and P. Michel for helpful discussions.

## References

  * Najmudin  _et al._ [2001] Z. Najmudin, M. Tatarakis, A. Pukhov, E. L. Clark, R. J. Clarke, A. E. Dangor, J. Faure, V. Malka, D. Neely, M. I. K. Santala, and K. Krushelnick, Measurements of the inverse faraday effect from relativistic laser interactions with an underdense plasma, Phys. Rev. Lett. 87, 215004 (2001).
  * Pitaevkii [1967] L. P. Pitaevkii, Electric forces in a transparent dispersive medium, Sov. Phys. JETP. 12, 1008 (1967).
  * Haines [2001] M. G. Haines, Generation of an axial magnetic field from photon spin, Phys. Rev. Lett. 87, 135005 (2001).
  * Naseri  _et al._ [2010] N. Naseri, V. Y. Bychenkov, and W. Rozmus, Axial magnetic field generation by intense circularly polarized laser pulses in underdense plasmas, Physics of Plasmas 17, 083109 (2010).
  * Frolov [2004] A. A. Frolov, Excitation of magnetic fields by a circularly polarized laser pulse in a plasma channel, Plasma Physics Reports 30, 698 (2004).
  * Bychenkov and Tikhonchuk [1996] V. Y. Bychenkov and V. Tikhonchuk, Magnetic field generation by short ultraintense laser pulse in underdense plasmas, Laser and Particle Beams 14, 55–62 (1996).
  * Sheng and Meyer-ter Vehn [1996] Z. M. Sheng and J. Meyer-ter Vehn, Inverse faraday effect and propagation of circularly polarized intense laser beams in plasmas, Phys. Rev. E 54, 1833 (1996).
  * Lecz  _et al._ [2016] Z. Lecz, I. Konoplev, A. Seryi, and A. Andreev, Gigagauss solenoidal magnetic field inside bubbles excited in under-dense plasma, Scientific Reports 6, 36139 (2016).
  * Deschamps  _et al._ [1970] J. Deschamps, M. Fitaire, and M. Lagoutte, Inverse faraday effect in a plasma, Phys. Rev. Lett. 25, 1330 (1970).
  * Horovitz  _et al._ [1997] Y. Horovitz, S. Eliezer, A. Ludmirsky, Z. Henis, E. Moshe, R. Shpitalnik, and B. Arad, Measurements of inverse faraday effect and absorption of circularly polarized laser light in plasmas, Phys. Rev. Lett. 78, 1707 (1997).
  * Ali  _et al._ [2010] S. Ali, J. R. Davies, and J. T. Mendonca, Inverse faraday effect with linearly polarized laser pulses, Phys. Rev. Lett. 105, 035001 (2010).
  * Nuter  _et al._ [2018] R. Nuter, P. Korneev, I. Thiele, and V. Tikhonchuk, Plasma solenoid driven by a laser beam carrying an orbital angular momentum, Phys. Rev. E 98, 033211 (2018).
  * Nuter  _et al._ [2020] R. Nuter, P. Korneev, E. Dmitriev, I. Thiele, and V. T. Tikhonchuk, Gain of electron orbital angular momentum in a direct laser acceleration process, Phys. Rev. E 101, 053202 (2020).
  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Orbital angular momentum of light and the transformation of laguerre-gaussian laser modes, Phys. Rev. A 45, 8185 (1992).
  * Tikhonchuk  _et al._ [2020] V. Tikhonchuk, P. Korneev, E. Dmitriev, and R. Nuter, Numerical study of momentum and energy transfer in the interaction of a laser pulse carrying orbital angular momentum with electrons, High Energy Density Physics 37, 100863 (2020).
  * Longman  _et al._ [2020] A. Longman, C. Salgado, G. Zeraouli, J. I. Apinaniz, J. A. Pérez-Hernández, M. K. Eltahlawy, L. Volpe, and R. Fedosejevs, Off-axis spiral phase mirrors for generating high-intensity optical vortices, Opt. Lett. 45, 2187 (2020).
  * Wu and Wang [2017] D. Wu and J. W. Wang, Magetostatic amplifier with tunable maximum by twisted-light plasma interactions, Plasma Physics and Controlled Fusion 59, 095010 (2017).
  * Shi  _et al._ [2018] Y. Shi, J. Vieira, R. M. G. M. Trines, R. Bingham, B. F. Shen, and R. J. Kingham, Magnetic field generation in plasma waves driven by copropagating intense twisted lasers, Phys. Rev. Lett. 121, 145002 (2018).
  * Longman and Fedosejevs [2020] A. Longman and R. Fedosejevs, Optimal laguerre–gaussian modes for high-intensity optical vortices, J. Opt. Soc. Am. A 37, 841 (2020).
  * Barnett and Allen [1994] S. M. Barnett and L. Allen, Orbital angular momentum and nonparaxial light beams, Optics Communications 110, 670 (1994).
  * Longman and Fedosejevs [2017] A. Longman and R. Fedosejevs, Mode conversion efficiency to laguerre-gaussian oam modes using spiral phase optics, Opt. Express 25, 17382 (2017).
  * Mackenroth  _et al._ [2019] F. Mackenroth, A. R. Holkundkar, and H.-P. Schlenvoigt, Ultra-intense laser pulse characterization using ponderomotive electron scattering, New Journal of Physics 21, 123028 (2019).
  * Akou and Firouzjaei [2020] H. Akou and A. S. Firouzjaei, Direct electron bunch acceleration by laguerre–gauss laser pulse, Physics of Plasmas 27, 093102 (2020).
  * Shi  _et al._ [2021] Y. Shi, D. Blackman, D. Stutman, and A. Arefiev, Generation of ultrarelativistic monoenergetic electron bunches via a synergistic interaction of longitudinal electric and magnetic fields of a twisted laser, Phys. Rev. Lett. 126, 234801 (2021).
  * Miyazaki  _et al._ [2005] S. Miyazaki, S. Kawata, Q. Kong, K. Miyauchi, K. Sakai, S. Hasumi, R. Sonobe, and T. Kikuchi, Generation of a microelectron beam by an intense short pulse laser in the $tem(1,0)$ $tem(0,1)$ mode in vacuum, Journal of Physics D: Applied Physics 38, 1665 (2005).
  * Lu  _et al._ [2007] W. Lu, M. Tzoufras, C. Joshi, F. S. Tsung, W. B. Mori, J. Vieira, R. A. Fonseca, and L. O. Silva, Generating multi-gev electron bunches using single stage laser wakefield acceleration in a 3d nonlinear regime, Phys. Rev. ST Accel. Beams 10, 061301 (2007).
  * Decker  _et al._ [1996] C. D. Decker, W. B. Mori, K. Tzeng, and T. Katsouleas, The evolution of ultra‐intense, short‐pulse lasers in underdense plasmas, Physics of Plasmas 3, 2047 (1996).
  * Arber  _et al._ [2015] T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, Contemporary particle-in-cell approach to laser-plasma modelling, Plasma Physics and Controlled Fusion 57, 113001 (2015).
  * Freidberg and Pearlstein [1978] J. P. Freidberg and L. D. Pearlstein, Rotational instabilities in a theta pinch, The Physics of Fluids 21, 1207 (1978).
  * Jackson [1999] J. D. Jackson,  _Classical electrodynamics; 3rd ed._ (Wiley, New York, NY, 1999).