The axion was originally proposed by Wilczek1 and Weinberg2 as a pseudoscalar particle emerging from the spontaneously broken Peccei–Quinn symmetry3,4 and providing an elegant solution to the problem of why the strong interaction conserves charge—parity symmetry (the strong CP problem).5 Generic models have predicted an inverse relation between the axion mass m_a and the coupling g_aγγ between the axion and electromagnetic fields.6–8 More generally, axion-like particles (ALPs) unrelated to the Strong CP problem may naturally arise from various theoretical frameworks beyond the Standard Model.9–14 Axions and ALPs are also prominent candidates for the cold dark matter in the Universe.15–18

The interaction between the axion and photon fields has been the focus of experimental efforts for decades. Axion haloscope experiments aim to detect interactions of galactic axion dark matter and magnetic fields via microwave cavities or magnetic flux loops.19–22 Helioscope experiments are axion telescopes to convert astrophysically produced axions (e.g., solar axions) into x rays in magnetic fields.23,24 Many astrophysical observations can also be used to search for and constrain the properties of axions.25–28 Alternatively, axions can also be directly produced and searched for in laboratories, particularly with intense laser beams. Rotation of the polarization of a laser beam in a transverse magnetic field, known as vacuum birefringence (VB), can arise from nonlinear quantum electrodynamics (QED) in strong fields.29 Signals deviating from the standard QED prediction would be a signature of the axion field. PVLAS,30 BMV,31 BFRT,32 and Q&A33 are experiments along this direction. Another representative class of experiments are based on so-called “light shining through a wall.” A laser is shone through a strong magnetic field with a downstream wall to block all of the laser beam, and if axions are produced, they will be transmitted through the wall and have a finite probability of being reconverted into photons in another magnetic field on the other hand of the wall. Such a scheme is used in experiments such as ALPS and ALPS-II,34,35 GammeV,36 LIPSS,37 OSQAR,38 and BMV.39

Owing to rapid developments in high-power laser technology and the related progress in producing ultrahigh magnetic fields in a laser–plasma environment, enhancement of polarization rotation has recently been reconsidered.40–42 An initially linearly polarized laser may even gain a small ellipticity after propagating through a magnetic field. Other axion-related phenomena may also be observable in laser–plasma interactions.43–45 However, current studies of axion-coupled laser–plasma interactions are still limited to ideal theoretical calculations.

In research into the standard laser–plasma interaction, particle-in-cell (PIC) simulation is a powerful tool, widely used in studying laser wakefield acceleration,46–49 strong-field QED effects,50,51 etc. In such a simulation, the particle motion is calculated according to the Lorentz force, with modifications from some QED processes, and electromagnetic (EM) fields are numerically solved according to Maxwell’s equations and the electric current contributed by the particles. In this paper, for the first time, we incorporate the axion field and its coupling with EM fields into the PIC code EPOCH.52 The axion field equation and modified Maxwell’s equations are numerically solved. Since the axion coupling is extremely feeble, the axion modulation to EM fields is considered as a first-order perturbation. Meanwhile, different types of axion field boundary conditions are introduced to satisfy different simulation scenarios. Such modifications allow us to simulate the axion field produced in a complex laser–plasma interaction self-consistently. We benchmark the PIC simulation with axion generation from photons and conversion into photons in a constant magnetic field against theoretical expectations.

In natural units (ℏ = ɛ₀ = c = 1), the axion coupling with EM fields can be described by the Lagrangian density term Lint=−14gaγγϕFμνGμν=gaγγϕE⋅B, where g_aγγ is the coupling constant, F and G are the electromagnetic tensor and dual tensor, respectively, and ϕ, E, and B are the axion field, electric field, and magnetic field, respectively. The total Lagrangian density readsL=LEM+Lϕ+Lint,(1)LEM=14FμνFμν−Aμjeμ,(2)Lϕ=12∂μϕ∂μϕ−12ma2ϕ2,(3)where m_a is the axion mass, and LEM and Lϕ are the Lagrangian densities of the EM fields and the free axion field, respectively. After variation of the total Lagrangian density, one gets the modified Maxwell equations with axion coupling:∂t2ϕ−∇2ϕ+ma2ϕ=gaγγE⋅B,(4)∇⋅E=ρ−gaγγB⋅∇ϕ,(5)∇⋅B=0,(6)∇×E=−∂tB,(7)∇×B=∂tE+j+gaγγ[(∂tϕ)B−E×∇ϕ],(8)where ρ and j are the electric charge density and current density, respectively. We have not introduced an additional duality symmetry here, as has been done in some other works.43,53,54 From Eq. (4), it is clear that the generation of axions requires that E · B ≠ 0.

We have implemented the numerical solver of the modified Maxwell’s equations in the EPOCH code, which is an open-source PIC code that is widely used in the laser–plasma community.52,55–57 EPOCH makes calculations in SI units, and in SI units, we choose the axion field ϕ to have dimensions of frequency and the axion coupling constant g_aγγ to have dimensions of time. In this way, the modified Maxwell equations in SI units read∂2c2∂t2−∇2+ma2c2ℏ2ϕ=gaγγℏμ0E⋅B,(9)∇⋅E=ρε0−cgaγγB⋅∇ϕ,(10)∇⋅B=0,(11)∇×E=−∂tB,(12)∇×B=∂c2∂tE+μ0j+gaγγc[(∂tϕ)B−E×∇ϕ].(13)

First of all, one needs to keep in mind that the axion modulation of the EM fields might be so small that it would be overwhelmed by the floating point error. To solve this problem, we use a method of field perturbation separation (FPS) by expanding the variables according to the order of axion coupling: E = E₀ + E₁, B = B₀ + B₁, ρ = ρ₀ + ρ₁, j = j₀ + j₁, where E₀, B₀, ρ₀, and j₀ satisfy Maxwell’s equations without the axion coupling. The first-order perturbations of the EM fields then satisfy∇⋅E1=ρ1ε0−cgaγγB0⋅∇ϕ,(14)∇⋅B1=0,(15)∇×E1=−∂tB1,(16)∇×B1=∂c2∂tE1+μ0j1+gaγγc[(∂tϕ)B0−E0×∇ϕ].(17)

As an estimation, we consider a laser with normalized vector potential a₀ = |eE₀/(m_eω₀c)| = 100, where ω₀ = ck₀ = 2πc/λ₀ is the laser frequency, and e and m_e are the electron charge and mass, respectively. This corresponds to an electric field E₀ = 4 × 10¹⁴ V/m for an 800 nm laser. Considering the generation of axions with such a laser field in a static magnetic field as strong as the laser magnetic field within a wavelength, the generated axion field is approximately gaγγE0B02μ0ℏk02=5.7×109s−1 (for a typical g_aγγ = 2.7 × 10⁻¹³ℏ GeV⁻¹). With such an axion field, the axion modulation of the EM field can be estimated as E₁ ≈ g_aγγϕB₀c = 4 × 10⁻¹³ V/m.

Unfortunately, the value of E₀ is represented as a double-precision floating-point number in the program. In binary, such a number is E0=bin1.011⋯0︸52×248 V/m. To represent such a number, the computer would use 1 bit for the ± sign, 11 bits for the exponent, and 52 bits for the mantissa. Therefore, the minimal step for E₀ to increase is error (E₀) = 2⁻⁵² × 2⁴⁸ = 0.0625 V/m, which is called the floating point error of E₀. We can see that E₁ is even far less than the floating point error: E₁ ≪error (E₀). Therefore, it is necessary to treat the zeroth- and first-order EM fields independently according to Eqs. (14)–(17). Otherwise, E₁ would be overwhelmed by the floating point error of E₀ and one could never get the correct modulation from the axion fields.

Moreover, since the Eq. (14) can be derived from Eq. (17) and the continuity equation ∂_tρ₁ + ∇ · j₁ = 0 together, we only need to treat the current perturbation j₁ in the code. For this part, we consider the perturbations of the particle momentum and velocity: p = p₀ + p₁ and v = v₀ + v₁. To calculate the perturbed particle motion, we should first calculate the perturbed fields felt by the particle, which are given byFx=x0+x1=(F0+F1)x=x0+x1=∑ijkS(rijk−x0−x1)(F0,ijk+F1,ijk)≈∑ijkS(rijk−x0)F0,ijk+∑ijk∂S(rijk−x)∂xx=x0⋅x1F0,ijk+∑ijkS(rijk−x0)F1,ijk=Fnop+Fp,(18)where F represents a field component such as E_x, E_y, E_z, B_x, B_y or B_z, S(r_ijk − x) is the shape function of a particle located at position x, r_ijk is the space grid position of the corresponding field component, and F_nop = ∑_ijkS(r_ijk − x₀)F_0,ijk and F_p are the unperturbed and perturbed fields, respectively, felt by the particle. The perturbed position and momentum of the particle can be calculated according todx1dt=v1,(19)dp1dt=q(Ep+v0×Bp+v1×Bnop),(20)v1=∂v∂p⋅p1=cγ03[γ02u1−u0(u0⋅u1)],(21)where q and m are the charge and mass of the particle, respectively, and u = p/mc is the normalized momentum. The current perturbation can then be calculated from the particle shape function as follows:(j0+j1)ijk=qΔV[S(rijk−x0−x1)(v0+v1)]≈qΔVS(rijk−x0)v0+qΔV∂S(rijk−x)∂xx=x0⋅x1v0+qΔVS(rijk−x0)v1,(22)where ΔV is the macroparticle volume in the code. So far, we have given a full description of the perturbed EM fields in the code, which is necessary to correctly calculate the derivative effects of axions.

In the following, we describe axion discretization in both the time and space grids. Superscripts will be used to indicate the index of time steps, and subscripts to label the space index. It should be pointed out that a single 0 or 1 as a subscript will be used to indicate the order of the EM field. Although we use SI units in the code, the following equations used in the numerical model will still be shown in natural units for convenience. From the equations, we naturally set the axion fields at the half-time grids. Then the time differential is realized as∂tϕ→(Dtϕ)n=ϕn+12−ϕn−12Δt,(23)∂t2ϕ→(Dt2ϕ)n+12=ϕn+32−2ϕn+12+ϕn−12Δt2,(24)where Δt is the time step. In each time step of the code, we will first update the EM fields for half a time step and then update the axion field according to the following difference equations:E1n+12=E1n+12Δt∇×B1n−j1n−gaγγ[(Dtϕ)nB0n−E0n×∇ϕn],(25)B1n+12=B1n−Δt2∇×E1n+12,(26)ϕn+32−2ϕn+12+ϕn−12(Δt)2=(∇2ϕ−ma2ϕ+gaγγE0⋅B0)n+12.(27)The particles are then pushed by both zeroth- and first-order fields. After that, we update the EM fields for the left half time step:B1n+1=B1n+12−12Δt∇×E1n+12,(28)E1n+1=E1n+12+12Δt∇×B1n+1−j1n+1−gaγγ[(Dtϕ)n+12B0n+1−E0n+1×∇ϕn+1].(29)Here, ϕ is given by the average of the axion fields at two time steps:ϕn=12ϕn+12+ϕn−12.(30)

For the spatial discretization, the configuration is shown in Fig. 1. While the components of the EM fields are staggered to the grid boundaries in certain directions, we consider the axion field to be locate at the cell center in all directions. In this way, we can calculate the space-dependent terms in Eqs. (25) and (27) as follows:∂x2ϕ→(Dx2ϕ)i′j′k′=ϕ(i′+1)j′k′−2ϕi′j′k′+ϕ(i′−1)j′k′(Δx)2,i′=i+12,j′=j+12,k′=k+12.(31)

Moreover, to update the axion fields as Eq. (27), we need to calculate the source term E₀ · B₀ at the axion field locations. The actual values of the EM fields are given by averaging the fields at the nearest cells:(E0⋅B0)i′j′k′=Ē0xB̄0x+Ē0yB̄0y+Ē0zB̄0z,(32)Ē0x=12E0x,i′+12j′k′+E0x,i′−12j′k′,(33)B̄0x=14B0x,i′j′+12k′+12+B0x,i′j′+12k′−12+B0x,i′j′−12k′+12+B0x,i′j′−12k′−12.(34)The other components are spatially averaged similarly. The cross-term E₀ × ∇ϕ in Eq. (25) is also given by averaging. Here, we show the expression for one component (E0×∇ϕ)x as an example:(E0×∇ϕ)x→E0yDzϕ−E0zDyϕi′+12j′k′,(35)E0yDzϕi′+12j′k′=14E0y,i′j′+12k′+E0y,i′j′−12k′×ϕi′j′(k′+1)−ϕi′j′(k′−1)2Δz+14E0y,(i′+1)j′+12k′+E0y,(i′+1)j′−12k′×ϕ(i′+1)j′(k′+1)−ϕ(i′+1)j′(k′−1)2Δz.(36)The other components are treated similarly.

The solutions of partial differential equations are only complete with specific boundary conditions. For the axions, we implemented four kinds of boundary conditions in the code: reflection, periodic, outflow, and perfectly matched layers (PML). These are quite similar to the treatments of the EM field in the usual PIC codes, which are used for different scenarios.58–601.Reflection boundaries. For the reflection boundary condition, we set the axion field gradient to zero at the boundary:ϕ−12j′k′=ϕ12j′k′,….(37)2.Periodic boundaries. For the periodic boundary condition, we set the axion field at the out cell to be the same on the other hand:ϕ−12j′k′=ϕNx−12j′k′,…,(38)where N_x is the total grid number in the x direction.3.Outflow boundaries. For the outflow boundary condition, we assume a left-handed wave propagating along the x direction at the x_min boundary, or a right-handed wave at the x_max boundary. The treatments in the y and z directions are similar. Specifically, we solve the two equations at the boundary to obtain the axion fields ϕ(−12)j′k′n and ϕ(−32)j′k′n−1 out of the simulation box at each time step, which are needed to update the axion field in the box. The results are as follows (with the source term S = g_aγγE₀ · B₀):ΔtΔx+(Δt)2(Δx)2ϕ−32j′k′n−1=−2ϕ−12j′k′n−1+2ϕ−12j′k′n−2−(Δt)2S−12j′k′n−1−ma2ϕ−12j′k′n−1+ΔtΔxϕ12j′k′n−1−(Δt)2(Δx)2ϕ12j′k′n−1−2ϕ−12j′k′n−1,(39)ϕ−12j′k′n=ϕ−12j′k′n−2+ΔtΔxϕ12j′k′n−1−ϕ−32j′k′n−1.(40)4.PML boundaries. For the PML boundary condition, we assume the axion field to decay exponentially in the absorbing layer:ϕ=ϕiei(ωt−kx)−σt,(41)where σ is an artificial absorption constant. Thus, we have performed the transform ∂_t → ∂_t + σ in Eq. (4). Then, after discretization, the evolution equation for the axion field in the absorbing layer can be written as(1+σΔt)ϕi′j′k′n+1=(Δt)2Si′j′k′n−m2ϕi′j′k′n+(Dx2+Dy2+Dz2)ϕi′j′k′n−σ2ϕi′j′k′n+σΔtϕi′j′k′n−1+2ϕi′j′k′n−ϕi′j′k′n−1.(42)

After implementing the axion field and the coupling with EM fields in the code, we can self-consistently simulate the axion generation, propagation, and conversion to EM fields. We first check the different kinds of boundary conditions in the code.

To check the boundary conditions used in the simulation, we assume a left-handed axion wave ϕ=e−(η2+y2)/w2⁡sink0η in the simulation box, where η = x + ct is the moving coordinate, w = 5λ₀ is the duration and width of the test axion pulse, and λ₀ = 2π/k₀ = 800 nm is the axion wavelength. To realize such an axion wave, what we actually do is to set the axion fields ϕt=12Δt and ϕt=−12Δt at two time steps. In this case, we consider the axion mass as 1 meV, which is much less than the mass corresponding to the axion wavelength set in the simulation. Since the here we only check the axion boundary conditions, the coupling constant g_aγγ does not matter. Moreover, the simulation box is set to be [−20λ₀, 20λ₀] × [−20λ₀, 20λ₀], and the grid number is 400 × 400. For the PML boundary condition, we set the layer thickness to be N_PML = 20 grids in all directions. The maximum absorption constant is σ_maxΔt = 0.2, and the absorption constant grows in the absorbing layer as σ=(1−iPML/NPML)3σmax, where i_PML is the grid index from the boundary.

Figure 2 shows the results obtained with different boundary conditions. We can see that the axion field can be correctly reflected or propagating to the other side with the corresponding reflection or periodic boundary conditions. With outflow or PML boundary conditions, the axion field can transmit out of the boundary or be absorbed by the boundary layer with negligible reflectivity. The numerical reflectivities are both ∼1% with the last two boundary conditions.

Moreover, we can further simulate the processes of photon–axion and axion–photon conversion in a constant magnetic field. Since the axion field equation contains the source term g_aγγE · B, a laser propagating in a constant magnetic field parallel to its polarization will continuously generate axions, and vice versa. The conversion probabilities are36,61Pγ→ϕ=Pϕ→γ≈14(gaγγB0l)2sin12κl12κl2(natural units),(43)Pγ→ϕ=Pϕ→γ≈c4μ0ℏ(gaγγB0l)2sin12κl12κl2(SI units),(44)where κ = k₁ − k₀ reflects the momentum difference between the photon and axion, k₀ and k₁ are the wave vectors of the photon and axion, respectively, and l is the laser propagation distance inside the magnetic field.

To simulate the process of axion generation, we consider a p-polarized laser with wavelength λ₀ = 800 nm propagating in a magnetic field B=Bcŷ. The normalized vector potential of the laser is a₀ = |eE/(m_eω₀c)| = 3, and it carries Gaussian envelopes in both transverse and longitudinal directions in the two-dimensional simulations. The width is 5λ₀ and the FWHM duration is 5T₀ = 5λ₀/c. The magnetic field is set to be 423 T and uniform in space. The laser propagates along the x direction, and the simulation window also moves at the speed of light after the laser has been injected into it. Moreover, here we assume the axion mass to be 1 meV. With regard to the coupling constant, there are two primary models for axion–photon coupling: the Kim–Shifman–Vainshtein–Zakharov (KSVZ) model6,7 and the Dine–Fischler–Srednicki–Zhitnisky (DFSZ) model.8,62 Different models will give different coupling constants:ma=6.3eV106GeVfa,(45)gaγγ=cγ2απfa,(46)where α is the fine structure constant, f_a is the energy scale of the axion, and c_γ = −0.97 for the KSVZ model and 0.36 for the DFSZ model. Here, the DFSZ axion model is considered. The corresponding coupling constant is g_aγγ = 2.65 × 10⁻¹³ℏ GeV⁻¹. The KSVZ model only gives a different coupling constant. The comparison between the simulation and theoretical results would still be similar.

Figure 3 shows the result of the axion generation process. The distributions of the laser electric field and axion field at two instants are shown in Figs. 3(a) and 3(b). The laser pulse continuously generates axions, and the axion field increases during the propagation. Meanwhile, we can also see that the laser pulse and the generated axion pulse both gradually defocus owing to the finite pulse width.

We also quantitatively compare the conversion probability. The conversion probability in the simulation is given by the total energy ratio between the axion and laser pulses in the simulation box:Pγ→ϕ,sim=HϕHγ=∬Hϕdxdy∬Hγdxdy,(47)Hϕ=ℏ2c1c2|∂tϕ|2+|∇ϕ|2+ma2c2ℏ2ϕ2,(48)Hγ=12ε0E2+1μ0B2.(49)The results are shown in Figs. 3(c) and 3(d). We can see that the conversion probabilities from theory and simulation match each other quite well, and Pγ→ϕ varies linearly with the propagation distance.

In addition, we can also simulate the process of axion conversion into photons. In this case, there is only an axion pulse initially in the simulation box and propagating in the constant magnetic field. The axion pulse is the same as that in Sec. III A, and the magnetic field is the same as that in Sec. III B.

In this case, the conversion probability in the simulation is also given by the energy ratio:Pϕ→γ,sim=HγHϕ=∬Hγdxdy∬Hϕdxdy,(50)where the energy density Hγ of EM fields is given by the first-order perturbations E₁ and B₁. The results are similar to those of photon conversion into axions, as shown in Fig. 4. The laser pulse is continuously generated by the axion pulse, and both pulses gradually defocus during propagation. The conversion probability obtained from the simulation again matches well with the theory.

As described in Sec. II B, we have used the FPS method to solve the problem of floating point error when there is a significant strength difference between the original and derivative fields. The validity of this method can also be checked by examining the propagation of an extremely weak laser pulse in a magnetized plasma. We consider a laser pulse with λ₀ = 800 nm and a plasma with a density n₀ = 3 × 10⁻³n_c, where nc=ε0meω02/e2 is the critical plasma density. The external magnetic field is set to be B_c = m_eω_p/e to magnetize the plasma in the z direction, where ωp2=n0e2/ε0m0 is the plasma frequency. The laser pulse is sufficiently weak, and the laser magnetic field is 10⁻²⁰B_c. The simulation box is [−20λ₀, 20λ₀] × [−50λ₀, 50λ₀], and the space grid number is 12 000 × 100. Higher spatial accuracy is adopted in the x direction, to reduce the impact of numerical dispersion. Such a laser pulse would be overwhelmed by the background magnetic field if it were calculated directly. However, its propagation can be correctly calculated with the FPS method.

The refractive index η of a laser in a transversely magnetized plasma is given by63ηo2=1−ωp2ω2,(51)ηe2=1−ωp2ω2ω2−ωp2ω2−ωp2−ωc2,(52)where the subscripts o and e indicate the ordinary and extraordinary waves, respectively, ω_c = eB_c/m_e is the electron cyclotron frequency in the magnetic field, and the magnetic field is perpendicular to the laser polarization. For weak laser propagation, we mainly consider the extraordinary wave. The laser group velocity is given byvg=cηeηe2+ω2d(ηe2)d(ω2).(53)On the other hand, the laser group velocity in the simulation is given by its centroid trajectory, and it can be calculated asxlaser=∬xB1z2dxdy∬B1z2dxdy.(54)

The results of weak laser propagation are shown in Fig. 5. Figure 5(a) shows that a sufficiently weak laser pulse can stably propagate in the magnetized plasma, rather than being numerically overwhelmed by the background magnetic field. Meanwhile, as Fig. 5(b) shows, the laser group velocity from the simulation matches well with the theoretical results.

In addition, the birefringence in the magnetized plasma is an interesting phenomenon, which arises from the different permittivities of different laser polarizations. Our FPS method can also handle such an effect correctly. To simulate the birefringent effect, we consider a higher plasma density n₀ = 0.33n_c to enhance the difference in refractive index between the ordinary and extraordinary waves. The external magnetic field is still as strong as B_c = m_eω_p/e. The weak laser pulse is now circularly polarized, and the peak magnetic field is 210−20Bc. We make the laser obliquely incident to the magnetized plasma, with an incident angle θ = π/6. After the laser has been refracted into the plasma, it is decomposed into two components: one is s-polarized and the other is p-polarized, as shown in Fig. 5(c).

To quantitatively check our FPS method as applied to the birefringence, we consider the wave vectors of the ordinary and extraordinary waves in the plasma, which are given bykx,o=k0ηo2−sin2⁡θ,kx,e=k0ηe2−sin2⁡θ.(55)With our parameters, we obtain k_x,o/k₀ = 0.65 and k_x,e/k₀ = 0.32, which match well with the Fourier transform in Fig. 5(d). Thus, our FPS method can also successfully reproduce the birefringent effect. In fact, not only can the FPS method solve the EM modulation from the axion fields, but it can also correctly handle weak lasers in strong background fields.

In summary, we have proposed and benchmarked a field perturbation separation (FPS) method to model the axion field and its coupling with EM fields in a PIC code. The axion field equation and modified Maxwell’s equations have been considered, with the modulation of the EM fields by the axion field being treated as a first-order perturbation to handle the huge orders of magnitude difference between the two types of field. Different boundary conditions are realized in the code. The conversion processes of axions into photons and of photons into axions have been checked as benchmarks of the code. The propagation and birefringence of a weak laser in a strongly magnetized plasma have also been checked to illustrate the validity of our FPS method. The new method and the PIC code based upon it may help researchers develop novel axion generation and detection schemes via laser–plasma interactions and achieve a better understanding of the astrophysical processes in which axions may participate.

Acknowledgment. The computations in this paper were run on the π 2.0 cluster supported by the Center for High Performance Computing at Shanghai Jiao Tong University. This work was supported by the National Natural Science Foundation of China (Grant Nos. 12225505 and 11991074).