Strong magnetic fields are essential for research in high-energy-density (HED) science, astrophysics, and controllable nuclear fusion. Plasma properties can be influenced by the presence of magnetic fields of different scales. For instance, astrophysical plasmas with a temperature and density of 10 eV < T < 100 eV and 10¹⁷ cm⁻³ < n < 10¹⁹ cm⁻³, respectively,1–5 can be magnetized with a magnetic field strength of about 10 T. Relativistic magnetic reconnection can be observed in such plasmas when the magnetic field strength reaches 1 kT,6–10 since the Alfvén velocity then becomes comparable to the speed of light. By contrast, the HED plasmas found in the cores of stars and planets, and in the fuel of inertial confinement fusion, have much higher temperatures and densities, exceeding 1 keV and 10²¹ cm⁻³, respectively. For these plasmas, a magnetic field strength of around 100 T can benefit electron transport when the cyclotron frequency is equivalent to or bigger than the collisional frequency. Additionally, a magnetic field strength of about 1 kT can be used to control relativistic electron beams, improving high-energy ion acceleration.11,12 When the magnetic field strength exceeds 10 kT, the laser can propagate in a magnetized plasma with electron density above the critical density as a whistler mode.13–16 The extremely small electron Larmor radius in magnetic fields above 10 kT leads to characteristic quantum-mechanical phenomena.17–19 The energy density of a magnetic field is given by ɛ_B = B²/(2μ₀), which can be expressed as ɛ_B [J m⁻³] ≈ 4 × 10¹¹B² [kT]. When the magnetic field strength exceeds 1 kT, we can study the properties and dynamics of matter under extreme conditions where the magnetic field energy density is greater than 10¹¹ J m⁻³. This may bring new opportunities to HED science.

The generation of strong magnetic fields in plasmas is of great interest owing to its potential for a wide range of applications,20 including magnetically enhanced fast-ignition fusion,21 generation of collisionless shocks in magnetized plasmas,22,23 magnetically assisted ion acceleration,24–29 and magnetic field reconnection research.30 To investigate the impact of strong magnetic fields, various approaches have been proposed and developed to generate fields over 100 T in laboratory environments. These include the use of pulsed-power devices,31–33 self-generated magnetic fields,34–38 magnetic flux compression,39–42 and laser-driven coils.43–45 The corresponding measurement techniques for the magnetic field have also made significant strides in development.46,47 The magnetic flux compression technique has become widely used for generating high magnetic fields in laboratory experiments. One advantage of this technique is the use of nondestructive, single-turn coils driven by pulsed-power devices to generate magnetic fields of several tens of tesla for long durations ($>$1 μs) and large volumes ($>$1 mm³).31–33 By utilizing these magnetic fields as seed fields, a strong magnetic field over 100 T can be generated using the flux compression method. Alternatively, the scheme based on a laser-driven coil can utilize high-power lasers with an intensity of ∼10¹⁵ W cm⁻² to generate a magnetic field exceeding 100 T without requiring compression and a seed magnetic field. However, the underlying mechanism by which the laser-driven coil produces such a strong magnetic field is not fully understood, and its effectiveness remains a topic of debate. In a recent study, Peebles et al.48 conducted extensive experiments to assess the potential of laser-driven coils in generating strong magnetic fields. They used different types of laser-driven coils to generate magnetic fields and employed various diagnostic techniques such as B-dot probing, Faraday rotation, and proton radiography49 to measure the fields. Their conclusion was that laser-driven coils cannot create quasi-static kT-level fields as claimed.48

At the same time, high-intensity laser systems have become increasingly important in strong magnetic field generation by inducing strong electric currents. With the invention of the chirped pulse amplification (CPA)50 technique, laser pulses can be amplified to possess an ultra-high energy density. A laser pulse with high intensity is an excellent driver for strong electrical currents suitable for generation of extremely strong magnetic fields. During the interaction between a high-intensity laser and an overdense plasma, the self-generated magnetic field on the plasma surface can be close in strength to the oscillating laser magnetic field.51 Although such high-amplitude magnetic fields are not in the bulk and have a complex topology, they can have an effect on generation of hot electrons.52

Generation of strong magnetic fields in a large volume requires the use of laser beams with both high intensity and high energy. Assuming the energy conversion from the laser into magnetic fields is 0.5%, a laser energy as high as 80 J is required to generate a 1 kT field with a volume of (100 μm)³. If the volume of the laser beam at focus is the same as the volume of the magnetic fields and the energy conversion is around 1%, we can expect that the strength of the generated magnetic field is going to be around B ≈ 0.1B_L,53 where B_L is the peak amplitude of the laser magnetic field. There has been a significant increase in the number of laser facilities around the world that are capable of producing peak power at the PW level.54 Meanwhile, the beam energy of several systems can be in the multi-kJ range. Such multi-kJ PW-class laser systems as LFEX,55 NIF ARC,56 and Petal,57 composed of multiple linearly polarized (LP) beamlets, allow experimental investigation of magnetic field generation in bulk plasmas within the relativistic regime by delivering the highest energy within picoseconds. The multibeamlet configuration is not only an essential feature of the laser system design, but also the key to advanced laser–plasma interaction regimes.12,58 The SG-II UP facility,59 which includes a ps PW laser and tens of kJ ns lasers, is now being upgraded to a larger-scale laser physics platform with multiple kJ-class ps laser beams and hundreds of kJ ns lasers. With the increased number of laser beams, the upcoming extended platform will provide the multibeamlet experimental capacity of kJ-class ps PW lasers.

One particular mechanism for producing self-generated magnetic fields is through the inverse Faraday effect (IFE), where an axial magnetic field is spontaneously generated when the laser transfers angular momentum (AM) to plasma electrons. There exist various methods for generating axial magnetic fields through AM transfer. In the early days, the IFE was mostly observed when circularly polarized (CP) radiation propagated through an unmagnetized plasma, resulting in the generation of a quasi-static axial magnetic field,60 but the effect is not exclusive to a CP beam carrying spin angular momentum (SAM) and is also applicable to a vortex beam carrying orbital angular momentum (OAM). It is now well known that helical wavefronts can be represented in a basis set of orthogonal Laguerre–Gaussian (LG) modes and that each LG mode is associated with a well-defined state of photon OAM.61 Magnetic field generation using LG beams has been investigated theoretically and numerically.62–64

While the IFE has been widely studied, the optical technology and elements required for producing CP or LG laser beams from conventional LP beams at high power can be costly, complex, and fragile. As of now, the generation of adequately strong and precisely controllable macroscopic fields at laser facilities employed in HED physics research54,59,65 remains a significant challenge. Recently, we have proposed a novel multibeam approach for AM transfer to a plasma, leading to subsequent strong magnetic field generation.66 The approach involves a spatial arrangement of four conventional laser beams, depicted in Fig. 1, that is inspired by the multibeam design of the PW-class kJ laser systems. The key aspect of our design is its ability to overcome the challenge faced by conventional LP Gaussian lasers, which is that they do not possess intrinsic AM and thus cannot achieve the IFE. By employing a specific spatial arrangement of laser beams, we enable the transfer of ensemble AM to the plasma, inducing strong rotating currents, which leads to efficient generation of a strong axial magnetic field. This new multibeam approach shows promising potential for enhancing magnetic field generation in laser-driven plasma systems.

Building on our previous work,66 this paper provides an in-depth analysis of the multibeam approach. We begin by examining the AM carried by several conventional laser beams in Sec. II to frame the study in the context of AM transfer. In Sec. III, we present the results of a 3D particle-in-cell (PIC) simulation for four laser beams with twisted pointing directions and a target containing a preplasma-like density ramp. A simple model of the axial magnetic field decay based on the expansion of hot electrons is provided. This section essentially reviews the key findings from Ref. 66, albeit for a different target, setting the stage for the subsequent discussion in Sec. IV. In Sec. IV, we examine various factors likely to impact field generation in actual experiments, including laser polarization, relative pulse delay, phase offset, laser pointing stability, and target configuration. Section V compares our scheme with other methods that use CP and LG laser beams to drive the magnetic field. Finally, Sec. VI summarizes the main results of this work.

It is well known that paraxial optical beams carry three distinct types of AM,67 namely, SAM for CP beams, intrinsic OAM (IOAM) for LG beams with helical wavefronts, and extrinsic OAM (EOAM) for normal beams propagating at a distance from the coordinate origin. The SAM is determined by the right- or left-hand CP. The IOAM is determined by the twist index of the LG beams. Both the SAM of a CP beam and the IOAM of a LG beam can be used for axial magnetic field generation via the IFE.60,62–64 We can understand the EOAM from the viewpoint of photons carrying linear momentum p at position r. The EOAM can be calculated as L^ext = r × p.

To take advantage of EOAM like SAM or IOAM in the IFE, we have proposed a setup that involves four laser beams with twisted pointing directions,66 which is shown schematically in Fig. 1. Each laser beam is represented by a colored cylindrical cone, with its direction determined by a wave vector k_i, where i is the index numbering the beam. The photon momentum in the ith beam is p_i = ℏk_i. For simplicity, consider two beams, the green and purple ones, with ${\mathit{k}}_{(1,2)}=(k_x,k_{\perp }^{(1,2)},0)$, intersecting the (y, z) plane at z_(1,2) = ±f₀ and y_(1,2) = 0, respectively, where f₀ is the beam offset. In this case, the axial AM of a given photon is given by [r × p]_x, where r is the position vector and p is the photon momentum. As a result, the total AM of the two beams is approximately$$L_x\approx -N\hslash (k_{\perp }^{(1)}-k_{\perp }^{(2)})f_0,$$(1)where N is the number of photons in one beam. If the two beams have the same tilt, then $k_{\perp }^{(1)}=k_{\perp }^{(2)}$ and, as a result, L_x ≈ 0. If the tilt of the second beam is opposite to that of the first, then $k_{\perp }^{(2)}=-k_{\perp }^{(1)}$. As a result, the total AM no longer vanishes, and it is roughly given by $L_x\approx -2N\hslash k_{\perp }^{(1)}{f}_0$. The total AM will double when a pair of such lasers with an offset in the y direction is added, showing that appropriately arranged beams can carry AM even though individual photons in each beam have no intrinsic AM.67 This observation bears similarities to studies involving γ-ray beams carrying OAM,68–70 where a population of photons with a twisted distribution of momentum p is generated. An alternative method is to calculate the AM of the electromagnetic field as$$\mathit{L}={\epsilon }_0\int \mathit{r}\times (\mathit{E}\times \mathit{B})d^{3}r,$$(2)where ɛ₀ is the dielectric permittivity and E and B are the electric and magnetic fields, respectively. Note that E × B matches the direction of the Poynting vector. In a given laser beam, the dominant component of the Poynting vector is directed along the axis of the beam. Therefore, four laser beams with a twist in the pointing direction (shown in Fig. 1) can possess net AM.

In this section, we present a self-consistent analysis performed using 3D PIC simulations with the open-source relativistic PIC code EPOCH.71 The key findings align with those first reported in Ref. 66. The overlap between parts of this section and Ref. 66 is intentional, serving to set the stage for the discussion in Sec. IV.

As explained in Sec. II, our goal is to take advantage of the multibeam configuration available at some state-of-the-art PW-class laser systems. We consider a laser system that provides four identical LP Gaussian laser beams that have no intrinsic AM. The peak intensity of each laser beam is I₀ = 8.0 × 10¹⁹ W cm⁻², which corresponds to a normalized electric field amplitude a₀ = 8.0. The normalization of the electric field is defined as a₀ = eE₀/m_eωc, where E₀ is the peak electric field amplitude, c is the speed of light, ω is the center frequency of the laser beam, and e and m_e are the electron charge and mass, respectively. The pulse duration τ_g = 600 fs, is the same for each beam (the temporal envelope of the electric field is Gaussian). The laser wavelength is λ = 1.053 μm and the beam waist radius is w₀ = 6.0 μm.

The orientation of the four laser beams is set according to Fig. 1(a), where the lasers are represented by colored conical cylinders. For each beam, the axis represents the propagation direction and the transverse size represents the beam radius. The lasers enter the simulation domain from the left boundary of the simulation box (x_e = −20 μm), which we refer to as the emitter plane. They interact with the plasma at the focal plane located at x_f = −5 μm. The transverse electric field of each beam in the (y, z) plane is set to be directed along the y axis. Figure 1(b) shows projections of all four beams onto the (y, z) plane. The beams are set up such that the intersection points of the beam axes with the (y, z) plane for a given longitudinal position x form the vertices of a square. This square formed by the intersection points rotates and shrinks as the beams propagate from the emitter plane to the target. To make it more evident that the intersection points get closer to each other as the beams approach the target, we have also plotted two circles that go through the intersection points in the emitter plane and the focal plane. The circle in the focal plane is visibly smaller. The radius of the circle in the focal plane is f₀. We call it the beam offset. The twist degree of the four laser beams is controlled by the azimuthal angle φ = 0.28π. By definition, there is no twist for φ = 0. The angle between the axis of every beam and the focal plane (in the plane formed by the beam axis and the x axis) is θ = arctan(S/D) = 0.27π, which is the polar angle characterizing the beam convergence, where D is the distance between the emitter plane and the focal plane and S is the transverse shift of the beam axes in the 2D projection plane.

We assume that the target is a fully ionized carbon plasma with an exponential longitudinal density profile to mimic a preplasma. The initial electron density is set to n(x) = 0.05n_e exp{(x [μm] + 5)/1.67}, where n_e = 50 n_c is the electron density of the foil, whose thickness is 3 μm behind the preplasma, and n_c = 1.0 × 10²¹ cm⁻³ is the critical density corresponding to a laser wavelength λ = 1.053 μm. Both electron and ion populations are initially cold (T = 0). The front surface of the foil is at x = 0 and the rear surface is at x = 3 μm. Using this profile, the plasma density ramps up from 3n_c to 50n_c over 5 μm, which increases the interaction volume between the laser beams and the plasma. The simulation box size is (40 μm)³, with grid cell sizes of (40 nm)³. We use four macroparticles per cell and have open boundaries throughout. Note that our resolution and the number of macroparticles per cell are comparable to what has been used in other studies of dense laser-irradiated targets and plasma-generated magnetic fields (see, e.g., Ref. 72). Table I gives detailed parameters of the 3D PIC simulation presented here. It must be pointed out that our setup differs from that used in Ref. 66, where we employed nanowires rather than a preplasma and the laser wavelength was set to 0.8 μm. We define the moment when the laser beams leave the simulation domain following their reflection off the target as t = 0 fs. In the described simulation, we observe generation and gradual evolution of an axial magnetic field. An illustration of the axial magnetic field B_x after the lasers have left the simulation box (t = 20 fs) is shown in Fig. 1(c). The blue, yellow, and red isosurfaces, from the outside to the inside, represent increasing magnetic field strength.

We observe generation of a strong axial magnetic field in the simulation when using beams with twisted pointing directions. Figure 2(a) displays the axial magnetic field for φ = 0.28π at t = 20 fs. The strength of the peak longitudinal magnetic field exceeds 10 kT. The volume occupied by the field stronger than 2 kT is ∼10⁴ μm³. The three volumetric isocontours indicate B_x/B₀ = −0.2, −0.55, and −0.8, where B₀ = 2πm_ec/|e|λ = 10.0 kT. To provide a clearer B_x profile, the values are temporally averaged over a 20 fs interval and spatially smoothed using a box with a stencil size of (0.4 μm)³. Owing to the approximately axisymmetric distribution of the magnetic field, for simplicity, we perform azimuthal averaging, which yields the magnetic field strength distribution in the (x, r) plane that is shown in Fig. 2(b). The magnetic field reaches its peak strength of 1.4B₀ on axis at x = −10 μm.

A strong axial magnetic field is present not only in front of the target, but also behind it. The magnetic field strength and its variation gradient are generally larger in front of the target. Although the magnetic field behind the target is weaker than the magnetic field at the front, its strength can still reach thousands of tesla. The magnetic field behind the target is generated by hot electrons carrying OAM that go through the target and exit at the rear side. The presence of the longitudinal field behind the target can be viewed as a confirmation that our scheme does indeed produce OAM-bearing electrons, since the lasers are unable to directly access this region.

To study the temporal evolution of the axial magnetic field, we perform longitudinal and azimuthal averaging over a cylindrical region of length L: $⟨B_x⟩(r)=(1/2\pi L){\int }_0^{2\pi }{\int }_0^L{B}_x(r,\varphi ,x)dxd\varphi$. We choose L = 5 μm according to the simulation result and show the distribution of ⟨B_x⟩(r, t) in Fig. 3(a). We can see from Fig. 3(b) that the magnetic field between t = −600 fs and t = −500 fs has the maximum growth rate Γ_in = 4.75 × 10⁻³ω, where ω is the laser frequency corresponding to the laser wavelength λ = 1.053 μm. As the hot electrons expand, the axial magnetic fields decay. We find that the characteristic decay rate is Γ_de = −2.4 × 10⁻³ω. The duration of the axial magnetic field is around τ_B ≈ 1 ps, which is of the same order of magnitude as the magnetic field generated by the IFE.64

Figure 4 shows the axial magnetic field distribution at three different times: panels (a)–(c) correspond to t = 20.0 fs, panels (d)–(f) correspond to t = 220.0 fs, and panels (g)–(i) correspond to t = 420.0 fs. Figures 4(a), 4(d), and 4(g) give the electron density normalized to the critical electron density n_c. The black contours denote n_e = n_c. The simulation confirms that the laser beams are reflected rather than being transmitted through the foil. Figures 4(b), 4(e), and 4(h) are transverse slices of the axial magnetic field (averaged over 20 fs output dumps and normalized to B₀ = 2πm_ec/|e|λ = 10.0 kT). The two black solid lines represent the front (x = 0 μm) and the back (x = 3 μm) of the foil, respectively. The front (x = −5 μm) of the exponentially modulated preplasma is represented by the black dashed lines. We find that the magnetic field region moves axially outward. Figures 4(c), 4(f), and 4(i) show radial lineouts of the axial magnetic field averaged azimuthally, temporally (over 20 fs), and longitudinally between the red and black dashed lines in Figs. 4(b), 4(e), and 4(h). The results show that the magnetic field expands radially outward, which leads to a reduction in its peak strength.

Although our main focus is on the axial magnetic field, we provide longitudinal slices of all three components of the generated magnetic field at t = 20 fs. In Fig. 5, B_x, B_y, and B_z are shown in the (x, y) plane (top row) and (x, z) plane (bottom row), respectively. These slices provide important reference values for the subsequent discussion of the dependence of the magnetic field on the twist angle. These results also confirm that the axial magnetic field is indeed very axisymmetric in our scheme. It is worth noting that in the region of peak axial magnetic field, the transverse components are generally an order of magnitude smaller in strength. Furthermore, these transverse fields reach their peak values away from the central core of the axial magnetic field. These transverse components are particularly weak near the axis. The described behavior can be clearly seen in Fig. 6, which provides the entire vector field structure and field lines in 3D at t = 20 fs. The videos in the Fig. 6 show the distribution from different viewpoints. We can conclude that the magnetic field is axial only in the central region.

To examine the efficiency of the magnetic field generation, we consider a box in front of the target with x ∈ (−15, −5) μm and y, z ∈ (−10, 10) μm. We find that the energy in the magnetic field is $U_B=\int (B_x^2/2{\mu }_0)dV\approx 1.6$ J. The total kinetic energy of plasma electrons within the same box is roughly ten times higher, U_e ≈ 20 J. The total incident energy in the four laser pulses is U_L ≈ 217 J. We thus conclude that the energy conversion efficiency from the lasers to hot electrons and from hot electrons to the plasma magnetic field are both around 10%. This overall conversion efficiency of ∼1% is similar to the conversion efficiency for the Biermann battery magnetic fields experimentally generated in laser–solid interactions.73 However, it is two orders of magnitude higher than the conversion efficiency for a laser-driven coil described in Ref. 74. This demonstrates the effectiveness of our multibeam approach in efficiently converting laser energy into a strong magnetic field.

We examined the kinetic electron energy E_k and the magnetic field energy U_Bx in a region in front of the target, defined by x ∈ (−20, −2) μm and y, z ∈ (−10, 10) μm. The time evolutions of these two energies are shown as red and blue lines in Fig. 7(a). The decay trends of both energies become nearly identical shortly after the laser pulses leave the simulation box (this occurs at t = 0 fs). This suggests that the observed decay of the magnetic field is primarily caused by the expansion of hot electrons, whose azimuthal current sustains the field. The expansion of hot electrons is influenced by several factors, including ion dynamics and potentially self-generated axial magnetic fields. A detailed study of the hot electron expansion will be addressed in future work.

To give a simple model of the expansion rate, we assume that the distribution of the axial magnetic field energy density has a Gaussian profile,$$u_B(x,y,z,t)=u_0\frac{V_0}{{\sigma }_x{\sigma }_r^2}\exp \left[-\frac{{(x-x_0)}^2}{2{\sigma }_x^2}-\frac{r^2}{2{\sigma }_r^2}\right],$$(3)$${\sigma }_x(t)={\sigma }_{x_0}+{\nu }_{x}t,{\sigma }_r(t)={\sigma }_{r_0}+{\nu }_{r}t,$$(4)with widths σ_x(t) and σ_r(t) along the different directions that vary over time to represent the expansion effects. The normalized factor $V_0/{\sigma }_x{\sigma }_r^2$ ensures that the total energy associated with the magnetic field, $U_B=\int u_B(x,y,z,t)dV=\int (B_x^2/2{\mu }_0)dV$, remains constant during the expansion. The peak of the axial magnetic field is assumed to be at x = x₀, r = 0. We determined the parameters in Eq. (3) through several steps. First, we used simulation data at t = 20 fs to calculate the energy of the axial magnetic field as a function of x and r, assuming axial symmetry. From the resulting profile, we found σ_x(t = 20 fs) = 3.5 μm, σ_r(t = 20 fs) = 4.5 μμm, and u₀V₀ = 4.4 mJ. We then scanned the parameters ν_x and ν_r to achieve the best fit to the blue line in Fig. 7(a), with the fitting result shown as a black dashed line with cross marks. The fitting parameters were determined to be ν_x = 2.7 × 10⁻² μm fs⁻¹ ∼ 0.09c and ν_r = 4.5 × 10⁻³ μm fs⁻¹ ∼ 0.015c. Consequently, we obtained σ_x [μm] = 3.0 + 0.027t [fs] and σ_r [μm] = 4.4 + 0.0045t [fs]. Our model shows that the expansion speeds differ significantly between the longitudinal and transverse directions, which can be partially attributed to the confinement by the generated axial magnetic fields.

In the proposed expansion model, the axial magnetic field, whose energy density is described by Eq. (3), is given by$$B(x,y,z,t)=\sqrt{2{\mu }_0\frac{u_0{V}_0}{{\sigma }_x{\sigma }_r^2}}\exp \left[-\frac{{(x-x_0)}^2}{4{\sigma }_x^2}-\frac{r^2}{4{\sigma }_r^2}\right].$$(5)After substituting all the parameters, we obtain the following explicit form:$$\begin{array}{ll}\hfill B_x(x,y,z,t)[\mathrm{k}\mathrm{T}]& =\frac{105}{(4.4+0.0045t[\mathrm{f}\mathrm{s}])\sqrt{3.0+0.027t[\mathrm{f}\mathrm{s}]}}\hfill \\ \hfill & \times \exp \left[-{\left(\frac{(x-x_0)[\mathrm{\mu }\mathrm{m}]}{6.0+0.054t[\mathrm{f}\mathrm{s}]}\right)}^2\right.\hfill \\ \hfill & -\left.{\left(\frac{r[\mathrm{\mu }\mathrm{m}]}{8.8+0.009t[\mathrm{f}\mathrm{s}]}\right)}^2\right].\hfill \end{array}$$(6)For the peak of the axial magnetic field $B_x^m$, the decay scaling obtained from our model is$$B_x^m[\mathrm{k}\mathrm{T}]=105{(4.4+0.0045t[\mathrm{f}\mathrm{s}])}^{-1}{(3.0+0.027t[\mathrm{f}\mathrm{s}])}^{-0.5}.$$(7)The result of Eq. (7) is represented by the black dashed line with cross marks in Fig. 7(b). The fitting model results in Figs. 7(a) and 7(b) are shown for t ranging from 0 fs to 1700 fs. Since the simulation data end at t = 920 fs, our model provides a longer prediction of the magnetic decay.

We begin by analyzing the azimuthal current density j_ϕ, which is believed to be responsible for the generation of the axial magnetic field. Figure 8(a) illustrates the distribution of the azimuthal current density j_ϕ (averaged over the azimuthal angle) for φ = 0.28π at t = 20 fs as a function of axial x and radial r positions. The current density is normalized to j₀ = −|e|cn_c = −4.8 × 10¹⁶ A m⁻² which corresponds to electrons with density n_e = n_c = 1.0 × 10²¹ cm⁻³ rotating at the speed of light. The establishment of the strong azimuthal current j_ϕ shown in Fig. 8(a) is crucial for the efficient generation of the strong axial magnetic field in our multibeam approach. The negative azimuthal current j_ϕ < 0 or j_ϕ/j₀ > 0, should produce a negative axial magnetic field (note that j₀ < 0). This is the same direction as shown in Fig. 2(a). 2D transverse slices of j_ϕ at x = −10 and −5 μm are presented in Figs. 8(b) and 8(c), respectively. As can be seen from Fig. 8(a), the azimuthal current density j_ϕ at x = −5 μm is stronger than the current density at x = −10 μm.

To estimate the maximum value of |B_x|, we make an order-of-magnitude assumption that j_ϕ is uniform within a cylindrical region of radius R and length 2h. After applying the Biot–Savart law,75 we obtain$$\begin{array}{ll}\hfill \max |B_x|& \approx \frac{{\mu }_0}{2}{\int }_0^R{\int }_{-h}^h|j_{\varphi }|\frac{r^2}{{(r^2+x^2)}^{3/2}}dxdr\hfill \\ \hfill & ={\mu }_0|j_{\varphi }|harsinh(R/h),\hfill \end{array}$$(8)where μ₀ = 1.26 × 10⁻⁶ H/m is the permeability of free space. From the information provided in Fig. 8(a), we can set R ≈ h ≈ 5 μm. Considering that the value of the azimuthal current density in Fig. 8(c) is |j_ϕ| ≈ 0.08|j₀|, where j₀ ≡ −|e|cn_c, we can estimate the maximum value of |B_x| as max |B_x| ≃ 20 kT. This result is close to the peak magnetic field strength ⟨B_x⟩ ≃ 1.8 B₀ = 18 kT observed in Fig. 3(b).

Assuming that the motion of the electrons is the main contributor to the current, we can calculate the effective azimuthal velocity as v_ϕ ≈ −j_ϕ/|e|n_e. Furthermore, the azimuthal current density can be used to estimate the OAM density of hot electrons. With the electron density obtained from the simulations, n_e ≈ 10²⁷ m⁻³, we find that the rotation velocity is about v_ϕ ≈ 0.1c, which implies a rapidly rotating plasma environment. One may wonder if the rotational effect is just a combination of four cross currents forming a shape similar to a square. We can estimate the expanding effect of the cross currents (if these currents are indeed present) using v_ϕ ≈ 0.1c. After 400 fs (from t = 20 fs to t = 420 fs), the transverse shift should be R_S = 0.1c × 400 fs ≈ 12 μm. Such a significant shift is not observed in our simulation, which leads us to conclude that we are indeed dealing with a rapidly rotating plasma rather than with four cross currents. We can write the OAM density of electrons as L_xe ≈ rγ_am_en_ev_ϕ, where γ_a is the relativistic gamma factor and n_e is the electron density. We next take into account that j_ϕ ≈ |e|n_ev_ϕ to obtain that L_xe ≈ rγ_am_ej_ϕ/|e|. With r = f₀ and ${\gamma }_a\approx {(1+a_0^2)}^{1/2}\approx 8$, the OAM density is approximately L_xe ≈ 2.4 kg m⁻¹ s⁻¹. Compared with the case where twisted lasers are employed to generate OAM,76 our scheme produces a rotating plasma environment with electron density and rotation velocity that are two orders of magnitude higher.

To study the rotational effect induced in the plasma and to quantify the AM transfer from the four laser beams with a twist angle of φ = 0.28π in our simulation, we analyze the densities of axial OAM in electrons (L_xe) and ions (L_xi). Both quantities are normalized to a reference value L₀ ≡ n_cm_ecw₀ ≈ 1.65 kg m⁻¹ s⁻¹, where L₀ represents electrons with density n_c rotating with azimuthal velocity v_ϕ = c at a radial position with r = w₀. The analysis is conducted at a time t = 20 fs. In Fig. 9, the densities of axial OAM for electrons (top row) and ions (bottom row) are shown at t = 20 fs. Figures 9(a) and 9(b) show angle-averaged L_xe and L_xi, respectively, as functions of x and r. Further insight can be provided by two-dimensional transverse slices of L_xe and L_xi at two different axial locations, x = −10 μm [Figs. 9(c) and 9(d)] and x = −5 μm [Figs. 9(e) and 9(f)]. We would like to point out that the OAM density in Fig. C1 from the Appendix of our previous work66 is artificially high owing to a data processing error.

As can be seen from Figs. 9(a) and 9(b), the axial OAM density distribution range of ions in the longitudinal direction is smaller than that of electrons. This behavior could stem from the challenges involved in altering the state of motion of ions. While the axial OAM density of electrons remains positive and gradually decreases away from the target in most of the region (x < −7.5 μm), another intriguing phenomenon is the occurrence of negative axial OAM density for electrons near the preplasma (x ∈ [−7.5, −5] μm), indicating a reversal of electron rotation at larger radial positions. Over time, the hot electrons undergo radial and longitudinal expansion away from the target, which is likely to contribute to the reduction in axial magnetic field strength. According to Fig. 9, we can approximate the distribution ratio of the OAM density between electrons and ions as η_ei = L_xe/L_xi ≈ 1%. We attribute this observation to the significant mass disparity between ions and electrons (where m_i/m_e ≈ 10⁴). Despite the higher OAM density of the ions, the velocity of the ions is still much less than that of the electrons. Since the current does not depend on particle masses, these simulation results confirm our earlier assumption that the dominant contribution to the azimuthal current component comes from electron motion.

To construct an analytical model for this new scheme, we begin by examining the AM density of the four laser pulses. This AM density, denoted by L_x = ɛ₀[r × (E × B)]_x, vanishes on the axis and reaches its maximum value at the beam offset. In light of the extensive simulation results we presented earlier, we now turn our attention to the process of AM absorption. We envision this absorption occurring within a cylindrical region of radius R and length 2Δh. The lasers, first and foremost, impart AM to electrons. Through the electron current, a strong axial magnetic field is generated. During the generation phase, the evolving magnetic field introduces an azimuthal electric field E_ϕ. This electric field, in turn, facilitates the transfer of absorbed OAM from electrons to ions. The electrons and ions are assumed to rotate rigidly in a shell with angular velocities ζ_e and ζ_i, respectively. Since the particles are moving away from the center, we expect this approach to work well only for estimating initial peak values. We introduce OAM to describe the rotational motion of electrons and ions, denoted by L_e = I_eζ_e and L_i = I_iζ_i, respectively. Here, I_e is the moment of inertia for electrons and is given by I_e = πR⁴Δhm_en_e. Meanwhile, I_i represents the moment of inertia for ions and is calculated as I_i = (m_i/Zm_e)I_e, where Z is the nuclear charge number of the ions. These OAM estimates capture the rotational dynamics of the charged particles. The overall evolution of the OAM of electrons and ions can be described by the following equations:77,78$$I_e\frac{d{\zeta }_e}{dt}=M_{\text{abs}}-M_E,I_i\frac{d{\zeta }_i}{dt}=M_E,$$(9)where M_abs denotes the torque resulting from OAM absorption. The term M_E denotes the torque due to the azimuthal electric field E_ϕ and can be expressed as M_E = ∫|e|E_ϕ(r)rn_ed³r ≈ |e|E_ϕ(R)I_e/m_eR. The rotational motion of electrons generates a current density j_ϕe, which can be approximated as j_ϕe ≃ −|e|n_eζ_eR. Furthermore, we can estimate the magnitude of E_ϕ(R) ≈ −R∂B_x/∂t by using Eq. (8). As a result, we can express M_E as $M_E\simeq (d{\zeta }_e/dt)I_e^{\prime }$, where $I_e^{\prime }\equiv {\omega }_p^{2}R\mathrm{\Delta }harsinh(R/\mathrm{\Delta }h)I_e/c^2$, with ${\omega }_p={(n_e{e}^2/m_e{\epsilon }_0)}^{1/2}$ being the plasma frequency. According to the parameters above, we obtain the result $I_e^{\prime }\gg I_e$. With these considerations in mind, we arrive at the expression ${\zeta }_e(t)={(I_e+I_e^{\prime })}^{-1}{\int }_0^t{M}_{\text{abs}}(t^{\prime })d{t}^{\prime }$. This equation illustrates that the rotation of electrons closely follows the temporal profile of M_abs(t), akin to an effect of effective inertia. Given that $I_e^{\prime }\gg I_e$ under the conditions of our study, we find that M_E is approximately equal to M_abs. Consequently, we obtain $L_i\simeq {\int }_0^t{M}_{\text{abs}}(t^{\prime })d{t}^{\prime }\simeq L_e{I}_e^{\prime }/I_e\gg L_e$. This result reveals that the total OAM of ions greatly exceeds that of electrons, consistent with the simulation findings presented in Fig. 9.

The absorption of electromagnetic AM is proportional to energy absorption in general. As in other studies of IFE, the transfer of AM from the four laser beams to the electrons can be assessed through the application of AM conservation principles.60,62,64,66 The number of absorbed photons is N_abs = U_abs/ω, where U_abs represents the absorbed energy, approximately U_abs ≃ η_absU_L. Here, η_abs is the absorption coefficient of laser intensity absorbed over the axial distance, assumed to conform to η_abs = η₀ exp(x/κ). The absorbed OAM carried by laser photons is subsequently transferred to both electrons and ions, with the fraction of the OAM carried by electrons denoted by η_ei ≡ L_xe/L_xi. Since the most of absorbed OAM is transferred to ions, we can obtain the OAM density of electrons as L_xe ≃ η_eiL_abs. Specifically, the axial OAM density of electrons can be approximated as follows:$$\begin{array}{cc}\hfill L_{xe}(x,r)& \approx 0.75\eta {\eta }_{ei}\frac{I_0{\tau }_g}{\kappa c}\sin (\theta )\sin (\phi ){\mathcal{F}}_{xr},\hfill \\ \hfill & {\mathcal{F}}_{xr}=r\exp \left[\frac{x}{\kappa }-\frac{2{(r-f_0)}^2}{w_0^2}\right],\hfill \end{array}$$(10)where I₀ denotes the peak intensity of the incident laser pulses and τ_g denotes their duration. For the sake of simplicity, we deduce from the simulation results that $\eta ={\int }_{x_e}^0{\eta }_{\text{abs}}dx\simeq \kappa {\eta }_0\approx 0.1$ (κ ≈ 3 μm). We utilize r = 6 μm and x = 0 μm to determine the peak OAM density of the electrons, L_xe ≈ 1.75 kg m⁻¹ s⁻¹. This outcome falls within a similar range to the peak OAM density (L_xe ≈ 1.69 kg m⁻¹ s⁻¹) presented in Fig. 9. Furthermore, it closely aligns with the result (L_xe ≈ 1.5 kg m⁻¹ s⁻¹) obtained through utilization of the azimuthal current density j_ϕ as demonstrated in Fig. 8, considering the roughness of the model.

To simplify the analysis, we have made some approximations for the absorption coefficient η_abs, assuming it to be independent of the laser peak intensity I₀ and incidence angles (the polar angle θ and azimuth angle φ). In reality, the absorption mechanism is more intricate.79–81 By utilizing the OAM absorbed by electrons L_xe and the associated azimuthal current density, we thus estimate the final axial magnetic field as$$\left|\frac{B_x}{B_0}\right|\propto \left|\frac{j_{\varphi }}{j_0}\right|\propto \eta {\eta }_{ei}\frac{a_0^2}{{\gamma }_a}\frac{c{\tau }_g}{\kappa }\frac{{\mathcal{F}}_{xr}}{r}\sin (\theta )\sin (\phi ).$$(11)In accordance with Eq. (11), the manipulation of the axial magnetic field is achievable by altering the sign of the twist angle φ, a proposition validated in Sec. IV.

The purpose of this section is to examine various aspects of our scheme that are likely to be important during experimental implementation at multi-kJ PW-class laser facilities such as the upcoming upgrade of SG-II. The comprehensive analysis presented in this section is based on a series of 3D PIC simulations and it addresses such factors as twist angle, polarization direction, phase, delay, and plasma front structure. We want to point out that we use a target with a preplasma and laser beams with wavelength of 1.053 μm, which differs from the use of nanowires and a laser wavelength of 0.8 μm in our previous work.66

To confirm the significance of the twist angle φ, we conducted two additional 3D PIC simulations. These simulations encompassed twist angles of φ = 0 and φ = −0.28π, the latter being opposite in twist direction to the original setup. These additional simulations were performed while keeping all other parameters the same as listed in Table I.

Figures 10 and 11 present 2D longitudinal slices of all three components of the time-averaged (averaged over 20 fs output dumps) magnetic field at t = 20 fs with twist angles φ = −0.28π and φ = 0.0π, respectively. When employing laser beams with the opposite twist, the axial magnetic field is reversed, as depicted in Fig. 10. On comparing Figs. 5 and 10, it is evident that the axially symmetric axial magnetic field B_x is reversed, while the transverse magnetic fields B_y and B_z still exhibit antisymmetric features. Moreover, the directions of B_y in the (x, y) plane and the directions of B_z in the (x, z) plane are reversed in Fig. 10 compared with Fig. 5. In Fig. 11, the results suggest that in the absence of a twist, i.e., φ = 0, the axial magnetic field is substantially diminished, while the effect on the transverse magnetic fields is relatively less pronounced owing to the existence of a longitudinal current even in the absence of a twist.

The generation of transverse magnetic fields occurs in all three scenarios owing to the presence of an axial current induced by the laser pulses. These transverse fields are particularly weak at small radii, underscoring the dominance of the axial magnetic field in governing the magnetic field within the central region. This observation confirms the recognition of a robust mechanism responsible for the generation of axial magnetic fields. It also suggests that the twist parameter can act as a knob for controlling the magnetic field within the central region, similar to the role of the topological charge l for LG beams in the IFE.

The polarization direction may play a role in the AM transfer process of our scheme that employs multiple LP Gaussian beams. Intuitively, it may seem that the direction of the transverse laser electric field is crucial for the generation of the azimuthal current. However, it is worth noting that the laser electric field is oscillating, while the current generated by this scheme is quasi-static. The quasi-static current and magnetic field evolve on a much longer time scale than the laser period. The electric field oscillations make it difficult for the laser beams with specific polarization directions to generate the quasi-static azimuthal current with a preferred direction. However, the polarization direction can affect the absorption of laser energy and AM. On the basis of the absorption characteristics of the plasma, different polarization directions result in different degrees of laser energy absorption within the plasma.82–85 This variation could potentially affect the distribution of hot electrons generated within the plasma.

To investigate the impact of the direction of laser polarization, we performed two additional 3D PIC simulations where the polarization direction was deliberately changed compared with that used in the original simulation presented in Sec. III. Figure 12 shows the direction of the transverse electric field in the (y, z) plane for each beam, where Fig. 12(a) shows the original polarization and Figs. 12(b) and 12(c) show the polarization in the two additional simulations. In the original simulation discussed in Sec. III, we set the polarization of all four laser beams such that the transverse electric field of each beam in the (y, z) plane was directed along the y axis, as shown in Fig. 12(a). We will refer to this simulation as Locked LP₀. In the first additional simulation, we rotated the polarization direction in every other beam by π/2. The resulting polarization is shown in Fig. 12(b). We refer to this simulation as Locked LP₁. In the second additional simulation, we introduced a random rotation ξ ∈ (0, π) into the original polarization direction for each laser beam. The polarization used in the simulation is shown in Fig. 12(c). We refer to this simulation as Random LP. The purpose of this simulation is to examine whether the multibeam scheme remains effective even when the linear polarization directions are randomly disordered. It is important to note that in this article we are primarily focused on the effects of the combination of four laser beams with linear polarization on the generation of axial magnetic fields. We do not consider other polarization states. This choice is mainly due to the design of the multibeam LP laser setup used in most high-power laser systems.

Figure 12(d) shows the time evolution of the spatially averaged magnetic field for the three different polarization setups: Locked LP₀ (black), Locked LP₁ (red), and Random LP (blue). It is evident that the blue and black curves almost entirely overlap, whereas the red curve is slightly lower than the other two. On the basis of these results, we infer that co-polarized beams (beams with the same polarization direction) are slightly more favorable for generating the axial magnetic field in our scheme. The result also suggests that even if the polarization directions are somewhat disordered, this will not significantly affect the effectiveness of our scheme. From the perspective of AM absorption, the axial AM carried by the four laser beams, as given by Eq. (1), depends solely on the longitudinal component of the wave vector. Therefore, as long as the twisted pointing directions of the four lasers remain unchanged, polarization changes should have a minimal effect on AM transfer and, consequently, on the magnetic field generation mechanism. This result is promising for experimental implementation, potentially eliminating the need for deliberate adjustments of initial linear polarization directions.

In high-power laser experiments, time delay can pose several challenges that affect the accuracy and reproducibility of experimental outcomes. Specifically, a time delay can cause the interactions between the individual laser beams and the target to occur at different times. This may potentially impact the interaction dynamics. Given the challenges posed by laser synchronization difficulties, it is valuable to examine the role of delays in our scheme. To address this, we examined two different delay scenarios in simulations. In the first scenario, we introduced random delays of up to 250 fs for each of the four laser beams, representing the short-delay case. In the second scenario, the delays of four beams were set to 0, 150, 600, and 1000 fs. The delay settings for the two scenarios are illustrated in Figs. 13(a) and 13(b), where the duration of each laser pulse is 600 fs.

We find that under the short-delay condition in the original simulation, the interaction between the four laser beams and the plasma leads to generation of hot electrons over a greater spatial extent, implying a higher energy transfer efficiency and much faster expansion of the hot electron population. Consequently, in the context of the short-delay scenario, various factors, including the strength and spatial range of the axial magnetic field, exhibit a magnitude that is an order of magnitude higher than those observed under the long-delay scenario.

Figure 13(c) presents the temporal evolution of the spatially averaged axial magnetic field for the two considered delay scenarios. The field for the case without delay, i.e., the no-delay curve, is also given for reference. We find that the differences between the short-delay case and the no-delay case are minimal. However, the long delay noticeably reduces the strength of the generated magnetic field compared with the case without the delay. The peak value reaches 40% of the peak value without the delay. Despite this reduction, the magnetic field still reaches the kT range, demonstrating that our scheme remains effective under sufficiently long delay conditions. This provides encouraging prospects for experimental implementation.

An important takeaway from the long-delay simulation is that a strong magnetic field can be generated even if all four beams do not temporally overlap. In our case, only two beams overlap, as seen in Fig. 13(b). However, as shown by Eq. (1) and the corresponding discussion in Sec. II, two lasers with opposite orientations are sufficient to provide axial angular momentum. In our simulation, there is a significant period during which at least two lasers interact with the plasma target simultaneously, ensuring the transfer of angular momentum and thereby supporting magnetic field generation.

In multibeam laser experiments, the phase relationship between different laser beams can produce different synthetic effects. For example, multiple laser beams with a locked phase difference can produce an interference pattern and then improve the efficiency of laser energy conversion into hot electrons.58 To investigate the effect of phase locking between four beams in our scheme, we introduced random initial phases ψ_rd ∈ (0, 2π) for each laser beam. Figure 14(a) shows the time evolution of the spatially averaged magnetic field for two different cases: locked phase and random phase. We find that the introduction of random phases does not cause any significant changes. We can conclude that the phase relationship in the simulation did not have a significant effect on the generation of the axial magnetic field. No phase control is required for the combination of multiple laser pulses in our scheme.

Laser pointing fluctuations are inherent in any laser system, and so it is important to understand their impact. These fluctuations can affect our scheme by altering the beam offset f₀, the azimuthal angle φ, and the beam convergence angle θ for each beam. Our approach is to treat fluctuations in the beam offset as independent of those in φ and θ, allowing us to assess their relative importance.

To examine the role of beam offset fluctuations, we introduced random displacements Δx and Δy for each beam. We changed the position of a given beam hitting the target to (x + Δx, y + Δy), where x and y are the locations prescribed by our setup. The range for the fluctuations was set to |Δx| ≤ 4 μm and |Δy| ≤ 4 μm, whereas the original offset for each beam was $f_0=\sqrt{x^2+y^2}=11$μm. Note that we kept the angles of incidence unchanged while varying the displacement. Figure 15 provides a comparison of the plasma magnetic field in the simulation with fluctuations and that in the simulation without fluctuations. Despite significant fluctuations (Δf₀/f₀ ≈ 50%) the variation in the spatially averaged magnetic field is only around 10%. The characteristic time scale at t > 0 remains unaffected, with the kT-scale field persisting over a ps.

To examine the role of angular fluctuations, we introduced random variations Δθ and Δφ for each beam. We adjusted the angles for a given beam to θ + Δθ and φ + Δφ, where θ = 0.27π and φ = 0.28π are the convergence and twist angles prescribed by our setup. The range for the fluctuations was set to |Δθ| ≤ 0.1π and |Δφ| ≤ 0.1π. Note that we kept the beam displacement the same as in the original simulation. Figure 15 presents a comparison of the plasma magnetic field from the simulation with fluctuations and the original result. The change in the spatially averaged magnetic field is about 50%. Such sensitivity is not unexpected, because we are varying the twist angle, which, as shown in Sec. IV A, is key to magnetic field generation. Despite the reduction in field strength, the characteristic time scale at t > 0 remains unaffected, with the kT-scale field again persisting over a ps.

Our analysis suggests that angular fluctuations may have a more significant impact than displacement fluctuations. Given that the amplitude of the angular fluctuations in our simulation is ∼18°, the obtained result demonstrates the robustness of our setup. In actual experiments, both types of fluctuations are likely to be present. To fully understand their combined impact, it is essential to consider the specifics of the laser architecture and the experimental setup.

After examining the impact of various laser beam parameters, we now shift our focus to the influence of the target shape, specifically the density profile at the laser-irradiated surface. In all simulations presented in this paper, we used a solid target with an exponential preplasma, with parameters detailed in Table I. This choice is motivated by the fact that high-intensity kJ-class laser systems inherently have a prepulse that inevitably causes target ablation and generates a density gradient at the front (laser-irradiated) surface. For comparison, in our previous work,66 we used a solid target with initially intact nanowires to enhance laser absorption. A sufficiently thick preplasma can achieve this effect as well.

To examine the impact of the density profile at the laser-irradiated surface, we performed two additional simulations: one with initially intact nanowires and one with an extended lower-density shelf. The setup for the simulation with nanowires was similar to that used in our previous work.66 The wires were 5 μm long and 0.4 μm thick. Their electron density was the same as that in the target (n_e = 50n_c). The spacing between the wires was 2 μm. All other parameters were the same as in Table I. In the simulation with the shelf, we replaced the exponential density profile at x < 0 with a 5 μm shelf of density 12.5n_c. The temporal evolution of the spatially-averaged longitudinal plasma magnetic field in these simulations is presented in Fig. 14(b), alongside the curve for the original simulation with the exponential density profile.

Although the general trend is similar for all three curves, there are several important conclusions that we can draw from Fig. 14(b). First, we note that the nanowires used in our previous work66 are not essential for the magnetic field generation. Furthermore, the exponential density profile and the shelf generate a stronger field, with an improvement of nearly 10% according to the simulation data. We also find that the target with the density shelf sustains the field for a longer period of time. At t = 900 fs, its field strength is almost twice that of a target with the initially intact nanowires. These changes are likely associated with energy absorption and hot electron generation (see, e.g., Refs. 86–88), but a dedicated study is required to clearly identify the mechanisms at play. If achieving a low-density shelf is the goal, one might consider preheating the target in a controlled fashion,88 rather than relying on a prepulse to achieve this.

In previous work on the IFE, CP Gaussian beams were among the first to be studied. For a CP laser beam, each photon carries a maximum SAM of σ_zℏ, where σ_z = +1 for right-hand and σ_z = −1 for left-hand polarization. Later, vortex beams such as LG beams, which are known to carry OAM, were also used to study the IFE. For LG beams, each photon carries an OAM of (σ_z + l)ℏ, where l is the twist index. It is valuable to compare the magnetic field generation in our multibeam scheme with the magnetic field generation via the IFE using CP or LG beams. According to the model used in previous works60,62,64 and assuming that the laser absorption rate and the electron density are constant in space and time, the following expression for the axial magnetic field can be obtained (see Ref. 64 for more details)$$\begin{array}{l}\hfill B_z=\frac{2.25{\eta }_{\text{abs}}{I}_0\tau }{e{n}_e\omega w_0^2}|{\psi }_l{|}^2\left[l\left(\frac{|l|w_0^2}{r^2}-2\right)\right.\\ \hfill \left.+{\sigma }_z\left(\frac{|l{|}^2{w}_0^2}{r^2}-4|l|+\frac{4r^2}{w_0^2}-2\right)\right].\end{array}$$(12)We have used a Gaussian temporal pulse shape, with τ being the temporal full width at half maximum (FWHM). In Eq. (12), n_e is the electron density, η_abs is the laser energy absorption fraction per unit length, I₀ is the peak intensity and w₀ is the beam waist. Using reasonable assumptions, Longman and Fedosejevs64 obtained the absorption rate as η_abs = n_e/2n_ccτ and then estimated the peak axial magnetic field strength of the |σ_z| = 1, |l| = 0 and |σ_z| = 0, |l| = 1 modes for τ ≳ 100 fs and n_e ≳ 0.01 n_c as $|B{|}_{\max }\approx (P[\mathrm{T}\mathrm{W}]{\lambda }^3[\mathrm{\mu }\mathrm{m}]/w_0^4[\mathrm{\mu }\mathrm{m}])[10\mathrm{k}\mathrm{T}]$, where $P=0.75\pi I_0{w}_0^2$ is the laser power and η_abs ≈ 0.2 mm⁻¹ is the value for the underdense plasma. The axial magnetic field |B|_max shows a fourth-power dependence on the laser beam waist w₀. Meanwhile, under our high-density plasma conditions, the absorption rate is expected to reach η_abs ≃ η₀ ≈ 32.3 mm⁻¹. Despite having a high absorption rate in our scenario, the total absorption may be reduced owing to a limited interaction length of the laser absorption process caused by the presence of the high-density plasma that is opaque to the laser beams.

To make a quantitative comparison, we ran two 3D PIC simulations: one with a single CP Gaussian beam and one with a single LP LG beam. In these simulations, we set the laser parameters such that the laser energy was equal to the laser energy in the four beams from Table I. Additionally, we chose the beam waist radius to be equal to the beam offset (w₀ = f₀) in Fig. 1(b). This choice was made to ensure that the transverse extension of the generated magnetic field matched that for the four-beam setup. Note that previous studies indicated that different laser modes with different values of l and σ_z produce radial magnetic field distributions with the intense region confined within 1.5w₀, as shown in Ref. 64. It is worth pointing out that the magnetic field has a similar radial distribution for both the CP Gaussian beam (|σ_z| = 1, |l| = 0) and the LP LG beam (|σ_z| = 0, |l| = 1). Inserting these parameters into the expression for |B|_max above, we can estimate the axial magnetic field strength to be about 0.4 kT. In Fig. 16, the logarithmic plots show the time evolution of the spatially averaged axial magnetic field for the three schemes mentioned here. The black curve represents the scenario with four LP Gaussian beams. The red and blue curves represent the cases of a single LP LG and single CP Gaussian beam, respectively. It can be seen that the latter two have magnetic field strengths about two orders of magnitude lower than the former. In our simulations, the single CP Gaussian beam performs less effectively, while our multibeam scheme shows distinct advantages under the same energy condition.

We have presented a comprehensive computational study of a setup in which a strong axial magnetic field is generated in an interaction of multiple conventional laser beams with a target. The magnetic field strength in the considered setup reaches 10 kT, with the strong field occupying tens of thousands of cubic micrometers. The field persists on a ps time scale. Although a CP Gaussian laser beam and an LG laser beam can carry intrinsic AM, native beams at current high-power laser systems are not CP or LG beams, but rather conventional LP beams. Moreover, the design of kJ-level PW-class laser systems is such that they are composed of multiple LP beams. This consideration motivated the development of our scheme using four laser beams with a twist in the pointing direction (as shown in Fig. 1) to achieve net OAM.

The key role of the AM carried by these laser beams in initiating the magnetic field generation process has been established. The simulations have examined various aspects, including the distribution of the axial magnetic field, azimuthal current, and the electron and ion OAM densities. Dependences on factors such as twist angle, polarization direction, delay, phase, and target structure have also been extensively investigated. Our study has validated the effectiveness of this approach under various conditions, confirming its robustness and practical feasibility. In particular, the twist angle of the laser pulse emerges as a critical driver for maintaining the azimuthal plasma current that maintains the orientation of the magnetic field. Furthermore, the PIC simulations and supporting theory indicate that the twist angle serves as a convenient control knob for adjusting the direction and magnitude of the axial magnetic field.

When compared with other schemes using CP or LG beams, the multibeam configuration has several advantages. In particular, our approach requires only an LP Gaussian beam configuration, making it suitable for advanced high-power, high-intensity multibeam laser systems such as the upcoming major upgrade SG-II UP. Despite the challenges associated with pointing directions, the results underscore the feasibility of achieving a strong and sustained axial magnetic field using thoughtfully designed multibeam setups. These studies contribute significantly to the understanding of laser–plasma interactions and expand the capabilities of high-power laser systems. The results may provide new opportunities to study kT-scale magnetic fields where the magnetic field energy density is greater than 10¹¹ J m⁻³, which is usually the baseline in HED science. Other studies of strong-magnetic-field physics can also be expected.

The supplementary material videos show 3D magnetic vector field and magnetic field line structure from differentviewpoints.

Acknowledgment. Y.S. acknowledges support by the National Natural Science Foundation of China (Grant No. 12322513) and USTC Research Funds of the Double First-Class Initiative, CAS Project for Young Scientists in Basic Research (Grant No. YSBR060). A.A. was supported by the Office of Fusion Energy Sciences under Award No. DE-SC0023423. Simulations were performed with EPOCH (developed under UK EPSRC Grant Nos. EP/G054950/1, EP/G056803/1, EP/G055165/1, and EP/M022463/1). The computational center of USTC and Hefei Advanced Computing Center are acknowledged for computational support.