Magnetic Nanoparticle with a Disordered Shell (Reproduce Hysteresis Loop from article)

[Victor Hugo Sosa Fierro]

Hello! I'm trying to reproduce the hysteresis loop described in this article. However, my calculated Hc​ value is much lower than expected (see attached hysteresis loop for reference).

The code I am using is provided below. I have experimented with adjusting the edgesmooth parameter and ext_InterExchange() function, but these changes have had minimal impact. The values for Aex, K_u, and Msat are taken directly from Table 2 of the article. Any suggestions would be greatly appreciated!

// Set mesh grid and geometry
SetGridSize(12,12,12)
SetCellSize(0.5e-9,0.5e-9,0.5e-9)
FixDt=1e-15
Temp = 0
// edgesmooth = 7

//Set geometry
d1 := 5.5e-9
d2 := 0.5e-9
core := Ellipsoid(d1,d1,d1)
shell := Ellipsoid(d1+d2,d1+d2,d1+d2).intersect(core.inverse())
setgeom(core.add(shell))
saveas(geom, "core-shell")

//Define two regions, default region is 0
defregion(1, core)
defregion(2, shell)

// Set magnetic parameters for core
Msat_core := 7.74e5
Aex_core := 1.5e-11
Ku_core := 8.9e4
alpha_core := 0.01

msat.setregion(1,Msat_core)
Aex.setregion(1,Aex_core)
Ku1.setregion(1,Ku_core)
anisU.setRegion(1, vector(1, 0, 0))
alpha.setregion(1, alpha_core)

// Set magnetic parameter for shell
order := 0.8
msat.setregion(2,Msat_core*order)
Aex.setregion(2,Aex_core*order)
alpha.setregion(2, 0.005)

//Set the interaction between the two regions
//ext_InterExchange(1, 2, Aex_core)


//set initial magnetization in regions or shapes
//default = randommag
m.setregion(1, uniform(1, 0.1, 0))
m.setregion(2, randommag())
saveas(m, "initial_magnetization")

//Set and external field in x direction rotated 1 degree to break the simmetry
angle := pi/180
Bmax  := 3.56
B_ext=vector(Bmax*cos(angle),Bmax*sin(angle),0)
tableadd(B_ext)
relax() //works all the time

//set up the hysteresis loop parameters
H_min  := -Bmax
H_max  :=  Bmax
H_step :=  0.01

//ramp the field down
for H:=H_max;H>=H_min;H-=H_step{
    B_ext=vector(H*cos(angle),H*sin(angle),0)
    minimize() //only works if you are close to the minimum
    tablesave()
}

//ramp the field up again
for H:=H_min;H<=H_max;H+=H_step{
    B_ext=vector(H*cos(angle),H*sin(angle),0)
    minimize()
    tablesave()
}