Prediction of Deterministic All-Optical Switching of Ferromagnetic Thin Film by Ultrafast Optothermal and Optomagnetic Couplings

All-optical switching (AOS) of magnetization induced by ultrafast laser pulses is fundamentally interesting and promises unprecedented speed for magnetic data storage that is three orders of magnitudes faster than the current techniques. For ferrimagnetic material, the AOS is attributed to magnetic circular dichroism and angular momentum transfer between sublattices. Recently, ferromagnetic material is demonstrated in AOS under multiple pulses. Since the magnetic field needed to flip the ferromagnetic magnetization within femtosecond timescale is unphysically high, some theories hypothesized that there exists a prolonged magnetic field beyond the pulse duration in the switching process. This is intuitively inconsistent with the phenomenological explanation based on the light-induced magnetic field arising from the inverse Faraday effect (IFE). Here, we numerically study the AOS process and provide new insights into the long-standing paradox of the duration of the induced magnetic field. We show that the prolonged magnetic field duration originates from the ultrafast optothermal and optomagnetic coupling. Moreover, we numerically studied both single- and multiple-pulse AOS under different coupling strength between spins and the thermal bath in the macroscopic Fockker-Planck and Landau-Lifshitz-Bloch model. This numerical model may provide a guide to find suitable ferromagnetic materials for AOS.


S4
. Fitting of B i by the analytical expression. If the inverse Faraday effect is large compared with the M field, B i can be approximated by an analytical expression of convolution of laser pulse and exponential decay. Figure S3 shows that the fitting is very good. Figure S4. Fitting of τ 0 by the empirical expression. The decay time constant can be fitted by an empirical expression depending on the material electrical conductivity σ and the laser beam diameter D. Figure S4 shows that the empirical expression is accurate. Figure S5. T e and B i evolution over time. The edge of the laser beam has B i field larger than noise level during T e above the Curie temperature. This favors the edge-dominated flipping. The center temperature T e is above the Curie temperature for too long time and result in random magnetization. Figure S6. Edge-dominated flip with random initial magnetization. Edge can be flipped by 20-µm laser beam pulse with random initial magnetization. The side length of the plotted area is 15µm.

Two-Temperature Model
The transient response of the electron and lattice temperatures under the sub-picosecond laser irradiation are calculated as below S1, S2 .
The electron system absorbs the optical heating and its temperature peaks in a few hundreds of femtoseconds. Simultaneously electron system transfers its thermal energy to lattice system via the electron-lattice interaction. The electron-lattice system approximately reaches their local equilibrium at about 1 picosecond. Eventually, the thermal energy diffuses into its neighbors.

Macroscopic Fockker-Planck and Landau-Lifshitz-Bloch (LLB) model
The macroscopic LLB model is used to describe the transient response of the magnetization as follow where γ is the gyromagnetic ratio. L 1 and L 2 are the longitudinal and transverse kinetic coefficients.
The λ charaterizes the coupling strength between the spin and the thermal bath. Below the Curie temperature, Gilbert damping parameters α 1,2 can be expressed as, where m e is the equilibrium spin polarization at finite temperature.
The B eff has the form of 2 2 2 2 ,0 1 1 , 2 where B i is the induced magnetic flux density, B ani is the anisotropy field, In supplementary equation (4), L 1 , L 2 , and H J depend on the electron temperature T e and they will reflect the phase change when the electron temperature is above the Curie temperature T C . The spin precession plays a minor role in the AOS process S4 therefore we neglect the first and third term of supplementary equation (4). The anisotropy field is set to be zero in our analysis.

S9
The classical approach in plasma science calculates inverse Faraday effect S5 as where α is the magneto-optical susceptibility. The inverse Faraday effect serves as a transient source term j IFE =×H IFE during the magnetic field calculation. In this treatment, the ∂E/∂t term is neglected which is an acceptable approximation in AOS process.

Approximate expression of the induced magnetic field B i
The induced B field only explicitly depends on laser and material parameters and has the approximate analytical form as  and 0 ( ) u t is the unit step function. B 0 is the background B field. τ 0 is the decay time constant and t 0 is the time of the laser peak. The treatment of convolution is accurate because Maxwell's equations are linear with respect to the loop current j in the absence of material nonlinearity or variations in M field. We approximate the relaxation function g(t) as an exponential decay. Although the realistic decay time constant in g(t) also varies with time, this approximate can still effectively capture the B i relaxation at the initial stage of the AOS (Fig. S3). The magnetic disturbances can be modeled by adding more terms to the analytical expression of the induced magnetic field.
The laser duration is represented by w. F IFE is related to the laser fluence I laser by For the 10-nm thick ferromagnetic material (σ = 9.18×10 6 S/m) used in this study, C, m, and n are fitted from the numerical results as 3.16, 0.917 and 0.847 respectively. This correlation is accurate to the level of 6% for a range of σ (2.30-91.8 10 6 S/m) and D (2-20 μm) which is shown in Fig. S4.

Scalar Model
The scalar model is derived to simplify and reduce the simulation time. The electron temperature is calculated by the two-temperature model. Only the z component of the M field is considered and the analytical expression of B i is used. Figure 6 in the main text is generated by solving supplementary equations (1), (2), (12) and (15).

Numerical implementation
The detailed numerical studies in this work are carried out using commercial multi-physics software COMSOL. We used Heat Transfer in Solids, Magnetic Fields modules to capture thermal and magnetic responses respectively. The in-house module for material magnetization is constructed by using the generic mathematics model for the macroscopic LLB equation. The axial symmetric geometry is modeled under the cylindrical coordinates with minimum grid size about 3 nm.
The virtual ferromagnetic material parameters for the numerical model are list in Table S1.