We first derive the equations used in the field-driven version of Brandt method, which were not given when $Q\left(I_{m}, f\right)$ of cylinder was first computed in [7].

The studied cylinder is placed along the z axis at $r \leqslant r_{8}$. It is characterized by the PL $E(J)$ relation expressed by Eq. (6). We calculate Q at different values of $I_{c}, I_{m}, n$, and f by executing the steps described in [8], [9], [10], [11], [14], [15]. Dividing the cylinder into N coaxial tubes of equal cross-sectional area $\pi r_{s}^{2} / N$, each being centered at

and having uniform current density and electrical field, $J_i$ and $E_i$, respectively, the computation is started by calculating an $N×N$ matrix of components

which is defined for converting $J_{j}(t)$ to the flux $\Phi_{i j}(t)$ produced by it within unit length between $r=r_i$ and a reference radius $r_a$ by

The field-driven version is formulated based on the Ampère law and the integral expression of the Faraday law, $\oint \mathbf{E} \cdot d \mathbf{l}=-d \Phi / d t$, with

where $E_a$ is a uniform applied electrical field, whose meaning will be discussed in Section 4.

The components of the reciprocal of the matrix, written as $M_{i j}^{E \text { rec }}$, are used for calculating from $J_{i}(t)$ at $t=k \Delta t(k=0,1,2, \cdots), J_{i}^{k}$, to that at $t=(k+1) \Delta t, J_{i}^{k+1}$, by using N simultaneous equations for $i, j=1,2, \cdots, N$,

Where

is obtained from the Ohm law with nonlinear resistivity defined from Eq. (6),

And

is the sinusoidal applied electrical field in the field-driven version. The time step is set $\Delta t=T / N_{t}$, where $T=1 / f$ is the period and $N_{t} \sim 10^{4}$ is the computed time step number in one period, computation is performed using the MATLAB software and started at $k=0$ with $J_{i}^{k}=0(i=1,2, \cdots, N)$.

If current is a sinusoidal function of $t=(k+1) \Delta t$ as

quantity $E_{a}^{k+1}$ corresponding to$I^{k+1}$ is adjusted back-and-forth to have the difference between the computed and the given values of $I^{k+1}$ being within $\pm I_{m} / 1000$ [8], [9], [10], [11].

Having obtained $J_{i}^{k}(t)$ for $i=1,2, \cdots, N$, the loss $Q_k$ per unit length in kth time step is calculated by

The final loss Q per cycle per unit length is calculated by

where the sum is taken for one period starting at a certain time $t_{s}=k_{s} \Delta t$ for a stabilized Q to be finally obtained.

For the same cylinder with the same divisions as in the field-driven version, the current-driven version is formulated based on the Ampère law and the differential expression of the Faraday law, $\nabla \times \mathbf{E}=-\partial \mathbf{B} / \partial t$, with a driving current $I(t)$. We start to calculate the components of an $N×N$ matrix to convert $J_{j}(t)(j=1,2, \cdots, N)$ into flux density $B_{i}(t)$ in the azimuthal ϕ direction at $r_{i+1 / 2}(i=1,2, \cdots, N-1)$, defined by Eq. (7) with i replaced by i+1/2, and $I(t)$:

The Faraday law with a current $I(t)$ is then written as

For the case of sinusoidal $I(t)$ with $ I^{k+1}$ expressed by Eq. (15), $ d I(t), d J_{j}(t), d E_{i}(t)$, and $ d r_{i}$ are written $ \Delta I^{k} \equiv I^{k+1}-I^{k}, \Delta J_{j}^{k} \equiv J_{j}^{k+1}-J_{j}^{k}, \Delta E_{i}^{k} \equiv E_{i+1}^{k}-E_{i}^{k}$, and $ \Delta r_{i} \equiv r_{i+1}-r_{i}$, respectively, so that Eq. (20) is written

The components of the reciprocal of the matrix, written as $ M_{i j}^{I \mathrm{rec}}$, are used for calculating from $ J_{i}^{k}$ to $ J_{i}^{k+1}$ at given $ \Delta I^{k} / \Delta t$ by using N simultaneous equations for $ i, j=1,2, \cdots, N$,

The final Q is calculated by using Eqs. (16), (17).

The field and current-driven versions are best used under the condition of a given $E_a(t)$ and $I(t)$, respectively. When the field-driven version is used for a given $ I(t), E_{a}(t)$ has to be adjusted to meet the given $I(t)$ to a certain accuracy, so that more computation time is required with less accurate $I(t)$. Lai and Gu proposed a way to have $E_a(t)$ calculated after each time step of $J_i(t)$ calculation [18]. In the present case, we have from Eq. (11) that

Where

from which $E_{a}^{k+1}$ is calculated as

In this version with a given $I(t)$, an extra step is inserted to calculate $E_a$ for each time step of $J_i$ calculation, so that computation time is increased compared with the field-driven version with a given $E_a(t)$.

We propose another version to extend the $N×N$ matrix for the field-driven version into an $(N+1) \times(N+1)$ matrix with the $E_{a}(t)$ and $J_i(t)$ simultaneously calculated for a given $I(t)$, so that computation time can be much reduced. The components of the new matrix for converting $d J_{j} / d t$ and $E_a$ to $E_i$ and İ are

The components of the reciprocal of this matrix, written as $M_{i j}^{I E \text { rec }}$, are used for calculating from $J_{i}^{k}$ to $J_{i}^{k+1}$ and the corresponding $E_{a}^{k}$ by

The studied rectangular thin strip is placed along the z axis, and its width w is along the x axis and thickness $d≪w$, along the y axis. Defining sheet current density $K=Jd$, the PL relation Eq. (6) is now written

where $K_c$ is sheet critical current density. Dividing width w into N equal elements, each centered at

and having uniform current density and electrical field, $K_i$ and $E_i$, respectively, the computation is started by calculating an $N×N$ matrix of components

which is defined for converting $ K_{j}(t)$ to the flux $ \Phi_{i j}(t)$ produced by it within unit length between radius $ r=\left|x_{i}-x_{j}\right| \text {, }$ or $ w /(2 \pi N)$ if $j=i$, and a reference radius $r_a$ by

Eq. (31) was previously written in [8], [10], [9], [11] following Brandt [12], [5], [13] with a quantity $r_a$ added at the present, which will be explained below.

All the following steps are the same as those for cylinder, with J replaced by K.

For the same strip with the same divisions as in the field-driven version, the current-driven version is formulated based on the Ampère law and the differential expression of the Faraday law, $ \nabla \times \mathbf{E}=-\partial \mathbf{B} / \partial t$, with a driven current I(t). We start to calculate the components of an $N×N$ matrix to convert $ K_{j}(t)(j=1,2, \cdots, N)$ into flux density $ B_{i}(t)$ in the -y direction at $ x_{i+1 / 2}(i=1,2, \cdots, N-1)$ and $I(t)$:

The Faraday law with a current $I(t)$ is then written as

For the case of sinusoidal $I(t)$ with $ I^{k+1}$ expressed by Eq. (15), $ d I(t), d K_{j}(t), d E_{i}(t)$, and $ d x_{i}$ are written $ \Delta I^{k} \equiv I^{k+1}-I^{k}, \Delta K_{j}^{k} \equiv K_{j}^{k+1}-K_{j}^{k}, \Delta E_{i}^{k} \equiv E_{i+1}^{k}-E_{i}^{k}$, and $ \Delta x_{i} \equiv x_{i+1}-x_{i}$, respectively, so that Eq. (35) is written

It is converted as

where $ E_{j}^{k}$ is expressed by Eqs. (12), (13) with J replaced by K. The final Q is calculated by using Eqs. (16), (17).

The same as the cylinder case with Eqs. (23), (24), (25), where $ E_{j}^{k}$ is defined by $ K_{j}^{k} \rho\left(K_{j}^{k}\right)$.

The same as the cylinder case with $ E_{j}^{k}$ defined by $ K_{j}^{k} \rho\left(K_{j}^{k}\right)$ and

The applied field $E_a$ was first written by Rhyner and Yazawa et al. as $ -\nabla \phi $, where ϕ is the applied scalar potential, in their equations [15], [14] and later by Brandt as $ -d \mathbf{A}_{a}^{E} / d t $, where $ \mathbf{A}_{a}^{E} $ is the vector potential corresponding to the applied electrical field [16]. In order to clarify its physical meaning, we imagine, following Norris [3], that the current I carried by the studied long cylinder of radius $r_s$ returns through a coaxial perfect conductor tube (with zero resistance) of inner radius $r_a$ and that a generator gives an electromotive force (emf), which corresponds to a uniform applied field $ E_{a}(t) $ at $ r \leqslant r_{a} $ in the z direction, to drive I. In this case, the current in the ith tube $ I_{i}(t)=J_{i}(t) \pi r_{s}^{2} / N $ will produce flux $ \Phi_{i i}(t)=\mu_{0} I_{i}(t) \ln \left(r_{a} / r_{i}\right) /(2 \pi) $ in the azimuth φ direction between $r=r_i$ and $r_a$ per unit length, according to the Ampère law, with an emf between $r=r_i$ and $r_a$ per unit length, $ E_{a}(t)-E_{i}(t) $, induced by $ d \Phi_{i i}(t) / d t $, according to the Faraday law. Thus, we have $ E_{i}(t)-E_{a}(t) $ to be produced by $ \left(\mu_{0} r_{s}^{2} / 2 N\right) \ln \left(r_{i} / r_{\alpha}\right) d J_{i}(t) / d t, $, which corresponds to the term in Eq. (8) for $i=j$. For $i<j$ and $>j$, the relevant flux $ \Phi_{i j} $ produced by $J_j(t)$ occurs between $r=r_j$ and $r_a$ and between $r=r_i$ and $r_a$, respectively, which validates Eq. (8) for all values of i and j.

With the imagined current return tube, which will be named as the reference tube for simplicity, electromagnetic fields are confined within $ r \leqslant r_{a} $. Since the flux produced by current in the cylinder and so the emf induced by its change between any value of r and $r_a$ decrease continuously with decreasing $ r_{a}-\tau $, the value of field $ E(r, t) $ must approach $E_a(t)$ with increasing r to ra. Therefore, $ E_{a}(t) $ appearing in the field-driven version of Brandt method for the cylinder case is actually the boundary condition $ E\left(r=r_{a, t}, t\right) $ for the inner region of the reference tube, where electromagnetic fields are computed.

The Q of a cylinder of $r_s=0.5$ mm, $I_c=60$ A, and n=30 driven by sinusoidal $ E_{a}(t)=-E_{a m} \cos 2 \pi f t(f=5 \mathrm{~Hz}) $ at the same amplitude $ I_{m}=I_{c} $ of $I(t)$ has been computed with N=50 and different values of $r_a$. The resulting $ E_{a m} $ and the resulting Q normalized by Eq. (3) as q are shown as functions of $ r_{a} / r_{s} $ in Fig. 1(a) and (b), respectively, from which we see $ E_{a m} $ to increase linearly with increasing $ \lg \left(r_{a} / r_{s}\right) $ and q to decrease from 1.156 to 1.099 with increasing $ r_{a} / r_{s} $ from 1 to above 1000.

In terms of ac impedance, current $I(t)$ in the cylinder may be regarded as the field $E_a(t)$ applied on it divided by an impedance per unit length of the cylinder, and the impedance may be considered as the sum of a resistive part corresponding to the ac loss and an inductive part corresponding to the emf induced by $I(t)$, or as the sum of a current dependent impedance of the cylinder itself corresponding to the flux at $r⩽r_s$ and a constant reactance of the air inductance corresponding to the flux between $r=r_s$ and $r_a$.

From Eqs. (8), (9) and considering Eq. (15) to be an approximate $I(t)$, which at sinusoidal $E_a(t)$ is not sinusoidal in general, the inductive part of $E_a(t)$ compensating the emf induced by $ I(t), E_{a}^{\prime}(t)$, is calculated as

so that we have $ E_{a m}^{\prime}=0.003472$ V/m and 0 for $r_a/r_s$ =10000 and 1, respectively, calculated using $f=5$ Hz and $I_m=60$ A, and a straight line of $ E_{a m}^{\prime}$ vs $ \lg \left(r_{a} / r_{s}\right)$ to be parallel to the computed $ E_{a m}$ vs $ \lg \left(r_{a} / r_{s}\right)$ curve in Fig. 1(a).

When $r_a/r_s$ is very large, the total impedance is dominated by the constant reactance, so that sinusoidal $E_a(t)$ results in sinusoidal $I(t)$. Thus, the results in Fig. 1(b) mean that at a given $i_m=1,q$ is the largest when $ E\left(r=r_{s}, t\right)$ is sinusoidal and the smallest when $I(t)$ is sinusoidal.

Our computations show that at the same value of $ i_{m}, q$ for sinusoidal $E_a(t)$ to be larger than for sinusoidal $I(t)$ is a general rule, and at a fixed frequency, significant difference between both occurs at large n and $I_m$ when the cylinder is fully penetrated [27]. This is why results for n=30 and $i_m=1$ are chosen to be shown. Owing to the non-linear resistivity ρ proportional to $J^{n-1}$, as expressed by Eq. (13), ac loss is mainly contributed within the time intervals around $ I= \pm I_{m}$, where the maxima of $ \left|E\left(r=r_{s}, t\right)\right|$ are very large.

When $I(t)$ is sinusoidal, there are narrow $ \left|E\left(r=r_{s}, t\right)\right|$ peaks occurring at $ I(t) \approx \pm I_{m}$, whereas when sinusoidal $I(t)$ is replaced by sinusoidal $ E\left(r=r_{s}, t\right)$, the $ \left|E\left(r=r_{s}, t\right)\right|$ peaks are replaced by broad sinusoidal maxima, so that q is larger for the latter than the former, especially for n and $I_m$ to be large so that $ |E|$ increases quickly with $ |I|$ following the PL relation of $ E \propto I^{n}$ in the full penetration state.

Being the actual electrical field $ E\left(r=r_{a}, t\right), E_{a}(t)$ can be used for calculating Q from the energy flow density, i.e., Poynting’s vector $ \mathbf{S}(\mathbf{r}, t)=\mathbf{E}(\mathbf{r}, t) \times \mathbf{H}(\mathbf{r}, t)$, and noticing $ I(t)=2 \pi r_{a} H\left(r=r_{a}, t\right)$ by

where $ E_{a}^{k}$ is expressed by Eq. (14) and $ I^{k}=I(t=k \Delta t)=\sum_{i=1}^{N} I_{i}^{k}.$. Thus calculated q is plotted in Fig. 1(b) by dots, which are located inside the circles for q calculated using Eqs. (16), (17). We see both results to be identical.

As Eqs. (16), (17), Eq. (40) may be used for any waveforms of $E_a(t)$ and $I(t)$;$I^k$ is expressed by Eq. (15) when $I(t)$ is sinusoidal with $ E_{a}^{k}=E_{a}(t=k \Delta t)$. The agreement in q calculated in both ways may serve as a necessary condition for the correctness of numerical computation.

We obtain from Eqs. (8), (10) for a cylinder with reference radius ra,pra, and ra/p,p being a positive constant, that

from which we have the scaling relation

since the logarithmic factors in Eqs. (41), (43) are $ \ln \max \left(r_{i}, r_{j}\right)-\ln r_{a} \mp \ln p$, respectively.

The validity of Eq. (44) is general for any cylinder with a certain $E(J)$ at a given $I(t)$. For a cylinder of $r_s=0.5$ mm with PL $E(J)$ of $I_c=60$ A and n=30 at $ I(t)=I_{m} \sin 2 \pi f t$ of $ I_{m}=0.1 I_{c}$.1Ic and f=5 Hz, $ E_{a}(t / T)$ curves for $r_a/r_s$ =0.1,0.2,0.5,1,2,5, and 10 and $r_a/r_s$ =2,4, and 8 are shown in Figs. 2(a) and (b), respectively. Note that in (a) the curve for $r_a/r_s$ =1 is the average of the curves for $r_a/r_s$ =0.1 and 10, 0.2 and 5, and 0.5 and 2, and in (b) the curve for $r_a/r_s$ =4 is the average of the curves for $r_a/r_s$ =2 and 8.

As discussed in Section 4 A, $E_a(t)$ is physically the uniform emf applied to $r⩽r_a$ when $r_a⩾r_s$, but mathematically, it can be applied to $r=r_a>0$ and acting to anywhere as calculated from the same equations derived physically from the Ampère law and the integral expression of the Faraday law for the case of $r_a>r_s$. It is very interesting that when $r_a≪r_s$, $E_a(t)$ can be largely different from $E(t)$ inside the superconductor, with much larger amplitude and opposite waveform. Just like sinusoidal $I(t)$ to be realized by sinusoidal $E_a(t)$ driven calculation at $r_a/r_s ≫1$, it also can be realized at $r_a/r_s ≪1$, with a dominated unphysical negative inductive reactance.

It is necessary to recalculate the $q(n,i_m)$ at $I_cf=300$ AHz for cylinder given in [7], using the more advanced current-driven version. The new results are listed in Table 1 and shown in Fig. 3. In comparison with Table I in [7] computed by the field-driven version under the same conditions without data at n=40, we can see that the difference for most data is less than 1% except at very small $I_m$ or at large $I_m$ and n, it can be larger. Therefore, even if the advanced current-driven version is used, more careful work has to be done to decrease the discretization error by changing N and $N_t$ and making extrapolation to $N,N_t→∞$ before more accurate $q(i_m,n)$ is obtained.

The $q(t_s)/q$ vs $t_s/T$ curves computed for sinusoidal $E_a(t)$ and sinusoidal $I(t)$ using the field-driven and current-driven versions, respectively, are plotted in Fig. 4, Fig. 5 for a superconductor cylinder of $r_s$=0.5 mm, $I_c$=60 A, and n=5 and 30 at f=5 Hz and different values of $i_m$.

In Fig. 4(a) for N=100,n=5,im=0.1, and $r_a/r_s$ =1,1.4,2,5, and 105 computed under sinusoidal $E_a(t)$, we see that with increasing $t_s/T,q(t_s)/q$ oscillates around 1 with gradually decreased amplitudes and each amplitude increases with increasing $r_a/r_s$. The first maximum occurring at $t_s/T$≈0.43 increases from 1.0025 to 1.0125 to 1.0165 for $r_a/r_s$ =1,1.4, and 10⁵. 10⁵ was chosen as the high end of $r_a/r_s$. As can be seen from Fig. 4(a), the sinusoidal $E_a(t)$ and sinusoidal $I(t)$ gets closer with higher $r_a/r_s$ ratio. In one way of understanding, $r_a$ can be considered the radius of the returning path for the closed circuit. High $r_a/r_s$ ratio increases the inductance of the circuit and reduces both the phase lag due to the finite flux flow resistance and the nonlinearity of the circuit due to the loss characteristic of the superconductor. In pure linear inductive circuit case, the sinusoidal $E_a(t)$ and sinusoidal $I(t)$ is actually the same, considering the π/2 phase difference introduced in Eqs. (14), (15).

In Fig. 5(a) for N=100,n=5, and $i_m$=0.1,0.167,0.333,0.5,0.667,0.833, and 1 computed under sinusoidal $I(t)$, we see that the first maximum increases from 1.0166 to 1.0184 with increasing $I_m$ from 0.1 to 0.5 and then decreases, and at $i_m=1,q(t_s)/q=1$ occurs suddenly at $t_s/T⩾0.25$. This is also true for N=100,n=30, and $i_m=1$ as seen from Fig. 3(b), where the first maximum increases from 1.0014 to 1.0038 with $I_m$ increasing from 0.1 to 0.833.

The situation corresponding to Fig. 1 for N=50,n=30,im=1, and $r_a/r_s$ =1,1.4,2,5,10,50,102,103, and 104 with two values of 105 and 10100 added, computed under sinusoidal $E_a(t)$ shown in Fig. 4(b) is more complex. At $r_a/r_s =1,$q\left(t_{s}\right) / q$ approaches about 1 at $t_s/T=0.25$, just as the case under sinusoidal $I(t)$. The first maximum increases from 1 to 1.011 with increasing $r_a/r_s$ from 1 to 10, and then decreases to 1 with further increasing $r_a/r_s$ to 10¹⁰⁰. The first minimum decreases from 1 to 0.9964 with increasing $r_a/r_s$ from 1 to 100, and then increases to about 1 with further increasing $r_a/r_s$ to 10¹⁰⁰.

It can be noticed that for n=5 and $i_m=0.1$, the waveform for $r_a/r_s$ =10⁵ computed under sinusoidal $E_a(t)$ agrees quite well with that under sinusoidal I(t), as shown in Fig. 4(a) by the solid olive curve and dotted curve, and that for n=30 and $i_m=1$, the waveform for $r_a/r_s$ =10¹⁰⁰ computed under sinusoidal $E_a(t)$ agrees in average with that under sinusoidal $I(t)$, as seen in Fig. 5(b) between the green dotted curve and the blue dotted line.

In order to compare $q\left(t_{s}\right) / q$ vs $t_s/T$ of sinusoidal $I(t)$ for n=30 with its CS counterpart for n→∞, we plot in Fig. 6(a) the low-$t_s/T$ part of the data whose high-$t_s/T$ part are plotted in Fig. (b), to show how $q\left(t_{s}\right) / q$ for each value of $I_m$ increases from a minimum value to 1 with increasing $t_s/T$ from 0 to 0.25.

Using Eq. (3), sinusoidal $I(t)$ expressed in Eq. (1) is normalized as $ i(t)=i_{m} \sin (2 \pi t / T)$.

According to the critical-state model, when I is initially applied from the zero-field cooled state, $J_c$ penetrates from the surface to a depth $ r_{p}=r_{s}[1-i(t)]^{1 / 2}$, which decreases from $ r_{p}=r_{s} \text { to } r_{m}=r_{s}\left(1-i_{m}\right)^{1 / 2}$ with increasing I from 0 to $ I_{m} \leqslant I_{c}$. When I is reversed to $-I_m$ or $I_m$ after reaching $I_m$ or $ -I_{m},-J_{c}$ or $J_c$ penetrates from the surface, with $r_p$ decreasing also from $r_s$ to $r_m$. The nonzero electrical field E is induced at r between $r_s$ and $r_p$ by the flux change during $r_p$ decrease, with the same direction as I so that a loss occurs. The loss is a function of $r_p$ with a double value in the reverse than the initial interval, where the local J change at $r_p$ are $2J_c$ or $J_c$, respectively. Thus, q in Eq. (2) for the normalized loss per cycle per unit length as a function of $I_m$ may be used as twice the normalized loss for $r_p$ to decrease from $r_s$ to $r_m$, and for a given $r_m$, it can be used for any $i^∗$ that is not larger than $i_m$. Therefore, we write the function defined by Eq. (2) as $q\left(i^{*}\right)$, write q in $q\left(t_{s}\right) / q$ as $q\left(i_{m}\right)$, and $q\left(t_{s}\right)$ as the sum of $q_{1}\left(t_{s}\right), q_{2}\left(t_{s}\right)$, and $q_{3}\left(t_{s}\right)$, defined at t from $t_s$ to $T / 4, T / 4$ to $3T/4$, and $3T/4$ to $T+t_{s}$, respectively. Detailed consideration for the loss in each time interval leads to

where $i_{s} \equiv i\left(t_{s}\right)$ is defined at $0 \leqslant t_{s} \leqslant T / 4$ for simplicity.

The CS $q\left(t_{s}\right) / q$ vs $t_s/T$ curves corresponding to those in Fig. 6(a) for n=30 are calculated using Eq. (45) and plotted in Fig. 6(b). We see that both are qualitatively similar with a remarkable difference as follows. The CS $q(0) / q$ decreases with increasing $i_m$ from 0.1 to 1 continuously, but $q(0) / q$ at n=30 decreases with increasing $i_m$ from 0.1 to 0.833 and then increases with further increasing $i_m$ to 1. In fact, there is a minimum $q(0) / q$=0.7862 occurring at im≈0.86. This anomaly corresponds to a sudden drop of the first $q\left(t_{s}\right) / q$ peak to 1 when $i_m$ increases from 0.833 to 1 in Fig. 5(b). It should be related to the difference in the q vs $i_m$ curves between finite and infinite n; the curve is smooth around$i_m$=1 for the former but cut at $i_m$=1 for the latter. It could be interesting to further study the detailed physical processes for this, but at the present, we would like to study the reason for the oscillations of $q\left(t_{s}\right) / q$ at $t_s/T$⩾0.25 and finite values of n.

The CS $q\left(t_{s}\right) / q$ at $t_s/T$<0.25 is calculated in the above, with $q\left(t_{s}\right) / q=1$ occurring at $t_s/T$⩾0.25 owing to repeated identical reverse current penetrations for each period. The PL J increases continuously with E at finite n, as expressed by Eq. (6), which is different from the CS case, where J can be $±J_c$ and 0 only. This results in the oscillating $q\left(t_{s}\right) / q$ vs $t_s/T$ curves at finite n shown in Fig. 5, which is explained as follows.

In the CS case, the loss $q_{2}(T / 4)$ or $q_{3}(T / 4)$ in the reverse interval from $t=T/4$ to $3T/4$ or from $3T/4$ to $5T/4$ is twice of $q_{1}(0)$ in the initial interval from $t=0$ to $T/4$, and the loss in the first half of the reverse interval is significantly less than that of the second half. Therefore, the average E, which is proportional to loss power, in the initial interval is smaller than that in the second half of the reverse interval. In the PL case with finite n, there is an effective critical current density $J_{c, \text { eff }}$ that increases with increasing the average E to penetrate from the surface [27]. Thus, in the initial interval with increasing I from 0 to $I_m$, a positive $J_{c, \text { eff }}$ penetrates to a certain $r_m$, written $r_{m1}$, and in the following reverse interval with decreasing I from $I_m$ to $-I_m$, a larger negative $J_{c, \text { eff }}$ will penetrate to a larger $r_m$, written $r_{m2}$, leaving$J_{c, \text { eff }}>0$ still to occur at $r_{m 1}<r<r_{m 2}$ with a positive bias current $ΔI$. This will make the PL $q_{2}(T / 4)$ larger than $q\left(i_{m}\right) / 2$, since penetrated negative $J_{c, \text { eff }}$ has to result in a total I decrease of $2 I_{m}+\Delta I$, instead of $2 I_{m}$ for the stabilized state. The situation for the next reverse interval is the opposite, with the PL $q_3(T/4)$ smaller than $q\left(i_{m}\right) / 2$, corresponding to the penetrated positive $J_{c, \text { eff }}$ to result in a total I increase of $ 2 I_{m}-\Delta I$. Such processes will repeat for the following reverse intervals, which makes $q\left(t_{s}\right) / q$ to oscillate with $t_s/T$ continuously. Since ΔI itself will decrease exponentially with an n and $ i_m$-dependent time constant, $L_{\text {in }} / R_{\text {in }}, L_{\text {in }}$ and $R_{in}$ being intrinsic inductance and differential resistance at r around the initial penetration front relevant to $ ΔI$, the $q\left(t_{s}\right) / q$ oscillation amplitudes will decrease with time in accordance, as seen in Fig. 5. General speaking, $ R_{in}$ increases with increasing n, so that $ q\left(t_{s}\right) / q$ stabilization is slower at n=5 than at n=30, as seen in Fig. 5(a) and (b). At each value of n,$ R_{in}$ increases quickly when $I_m$ approaches 1, with $q\left(t_{s}\right) / q=1$ to occur accurately at $t_s/T⩾0.25$, just like the situation in the CS case.

More complicated $q\left(t_{s}\right) / q$ oscillations under sinusoidal $E_a(t)$ shown in Fig. 4 may be attributed to the contribution of an nonzero dc component of non-sinusoidal $I(t), T^{-1} \int_{t_{s}}^{t_{s}+T} I(t) d t$, to $ΔI$ and the contribution of $r_a/r_s$ -dependent external inductance $L_{ext}$ to the $I(t)$ waveform. For fixed n and $i_m$, the $q\left(t_{s}\right) / q$ vs $t_s/T$ curve approaches that of sinusoidal $I(t)$ when $r_a/r_s$ is large enough.

All results above are calculated for a cylinder with a clear physical meaning of $E_a(t)$ as described in Section 4.1. Different from a cylinder, $r_a$ for a strip is the inner radius of the reference tube coaxial with each j-th element in order to have a uniform $E_a(t)$ within the strip as written in Eq. (10), where J is replaced by K. In fact, $E_a(t)$ is the electrical field applied to $x=x_{j} \pm r_{a}(j=1,2, \cdots, N)$, it becomes the uniform applied field in the strip via the integral expression of the Faraday law, with $E_{i}(t)-E_{a}(t)=\left(\mu_{0} w / 2 \pi N\right) \ln \left(\left|x_{i}-x_{j}\right| / r_{a}\right) d K_{j}(t) / d t$ when $i \neq j$, for example, and it may be regarded as the boundary field at $r=r_a$ measured from the strip axis only when $r_a≫w$ with the small shifts of all the reference tubes along the x axis within ±w/2 to be negligible, so that E is practically uniform on the inner surface of the same reference tube of inner radius $r_a$.

In spite of these complexities, the features for the results of cylinder described in Sections 4.2, 4.3 Scaling of, 4.5 Stabilization of, 4.6 Critical-state, 4.7 Oscillating are still valid for strip at least approximately if $r_s$ is replaced by $w/2$, since the same Eq. (10) is used for both cases. In particular, Eq. (44) as the scaling relation is valid exactly for a strip as derived from Eqs. (32), (10); and Eq. (40) can be used for a strip conditionally, with Q to be calculated accurately only when $E_a(t)$ is defined at so large $r_a$ that actual $E(t)$ at $r=r_a$ becomes practically uniform.

When Brandt proposed his method for a strip or rectangular bar [5], [13], $r_a$ was set equal to 1 in the basic expression $\ln \left|\mathbf{r}-\mathbf{r}^{\prime}\right| /(2 \pi)$ of the matrix similar to that in Eq. (31), which was followed by most authors to calculate transport and magnetic properties [14], [15], [8], [10], [9], [11], [28], [29], [30], [31], [18], although the logarithmic of a quantity with a unit of meter seemed to be meaningless.

For the Q computation of a rectangular bar under sinusoidal $I(t)$ using the field-driven version, a length $r_0$ was included in all components of the matrix by Yazawa et al., who stated that as the contour $r=r_0$ gives an equal vector potential line, the constant r0 should be large enough compared to the conductor dimension, e.g., $r_0=1$ m [15].

On the other hand, when Brandt developed his method to compute transport properties of rectangular bars, he included the bar length L to all the matrix components but not for magnetic case, and stated that some electrodynamic properties of long bars with applied current (e.g., their self induction) depend logarithmically on the bar length L, while for bars in a magnetic field usually the limit $L→∞$ may be taken [16].

Introducing $r_0$ or L in the components of the matrix to make the argument of the logarithmic function dimensionless, their meaning is actually unclear.

In the present work for a simpler case of cylinder, a length $r_a$, equivalent to $r_0$ and L, is defined as the radius at which $E(t)=E_{a}(t)$ occurs.

Commenting on the statement in [15], we may say that Q under the same sinusoidal $I(t)$ expressed by Eq. (1) is independent of the value of $r_0$, so that $r_0$ to be large enough is not necessary, although such a large $r_0$ can be regarded as the radius at which $E(t)= E_a(t)$ occurs quite accurately even if the superconductor is not a cylinder.

Commenting on the statement in [16], we may say that Q at a given $I_m$ may decrease with increasing L if sinusoidal $E_a(t)$ is applied to the bars, as described in Section 4.2, and that the process to approach stabilization after applying a sinusoidal $E_a(t)$ will be influenced by the value of L owing to the change of dc component of I(t), as described in Section 4.7.

The length L of the superconductor to be used as the radial distance to its axis seems to mean that if a long superconductor of length L is elongated by connecting two long normal conductors in its each end, then the flux on the midplane of the superconductor produced by the current may follow the rules determined by the $E \propto J^{n}$ relation up to a radial distance L, beyond which the contribution of normal conductors with $E \propto J$ will be dominant. However, this assumption is too much idealized and it would be meaningless if the actual current flowing circuit has a radial size less than L and it will lead to $r_{a} \rightarrow \infty$ in Eq. (31) for the infinitely long strip, which is not the case from our analysis in the present work.

Since $M_{i j}^{E}$ expressed by Eq. (31) with $r_a=1$ m taken arbitrarily may be used not only for transport but also for magnetic case, as shown in [8], [10], [9], [11], [28], [29], [30], [31], the condition $L→∞$ stated in [16] for magnetic case is questionable. The absence of L in the matrix for magnetic case in [16] is a consequence of the antisymmetry of current in the bar under magnetic field; $\ln \left(\left|x_{i}-x_{j}\right| / r_{a}\right)$ in the matrix for a strip can be replaced by $\ln \left(\left|x_{i}-x_{j}\right| /\left|x_{i}+x_{j}\right|\right)$ with $N/2$ elements under a magnetic field. This idea was later emphasized by Brandt in [32], where instead of mentioning $L→∞$ in the magnetic case, he stated the strip length L to drop out from the symmetric kernel (the matrix in the present work).

In the field-driven version proposed by Brandt, the equation of motion for the current density was deduced in terms of vector potential $\mathbf{A}$ defined by $\mathbf{B}=\nabla \times \mathbf{A}$ and $\nabla \cdot \mathbf{A}=0$ and written for an infinitely long superconductor of cross-section $S^{\prime}$ as

where $r_a$ is added to its previous expression in [5], [13], [33] to make the dimensions correct. $E_a(t)$ is a gauge convariant quantity relating to the value of $r_a$ and not measurable in general with a complex meaning described above.

This difficulty is avoided in the current-driven version where all physical quantities are measurable. In this sense, the current-driven version is intrinsically better than the field-driven version.

Both the current and field-driven versions have nothing to do with $-∇ϕ$ produced by distributed electrical charges whose field obeys the Gauss law, so that $E_a$ should not be expressed by $-∇ϕ$ and the $\mathbf{A}-\phi$ formulation named in the literature may be questionable [33], [34].

Following Brandt [12], [13], previous computations of Q for strips or rectangular bars were all conducted using the field-driven version with $r_a$ replaced by 1 [14], [15], [8], [10], [9], [11], and in this case, in order to meet Eq. (1), $E_a(t)$ has to be back-and-forth adjusted, so that the computation is very slow even if each $I^k$ value is adjusted to an accuracy of $I_m/1000$.

We introduce a reference radius $r_a$ in the present work, and found that a sinusoidal $E_a(t)$ may result in a sinusoidal $I(t)$ when $r_a/r_s$ or $r_a/w$ is sufficiently larger or smaller than 1. Thus, with enough large $r_a$, we need only to adjust $E_{a,m}$ to get $I_m$ accurately, so that the computation time may be much reduced.

Such adjustments of $E_a(t)$ is omitted by the current-field-driven versions A and B, by which the given $I(t)$ is exactly set with the corresponding $E_a(t)$ obtained accurately. The $E_a(t)$ itself is omitted by the current-driven version, by which q is computed the quickest at the exact given $I(t)$.

To compare the computation speed of q at sinusoidal $I(t)$ among the current-driven and current-field-driven versions, we have made a test as follows. For a cylinder of $r_s$=0.5 mm, with PL $E(J)$ of $I_c$=60 A and n=30, when N=100,N_t=2×10⁴, and $r_a=r_s$, the same q=3.556×10^-4 is calculated at $t_s/T$=2,$i_m$=0.09, and f=5 Hz by the current-field-driven versions A and B and the current driven version with computation time of 7.41, 1.80, and 1.76 s, respectively. Thus, the current-driven version is the quickest, the current-field-driven version B is slightly slower but it can get $E_a(t)$ simultaneously, and having the same applicability as version B, the use of version A should be excluded for its 4 times more computation time.

It should be noticed that the minimum $t_s/T$ to get a stabilized q under sinusoidal $I(t)$ decreases with increasing n and im. As seen from Fig. 5(a) and (b), $t_s/T$=2 and 0.7 are large enough to have $q\left(t_{s}\right) / q$ deviating from 1 within ±0.1% at n=5 and 30, respectively, for any value of $i_m$, and the deviation can be less than ±0.01% at $i_m=1$ and any value of n when $t_s/T$ is as small as 0.25.

We should mention that the current-driven version can also be used to solve magnetic problems like the field-driven version. For the studied current carrying strip under a perpendicular magnetic field $H_{u}(x, t)=B_{a}(x, t) / \mu_{0}$ applied in the y direction. Faraday’s law is written as

where $\Delta B_{\mathrm{ai}}$ is the external magnetic field at (x_i+x_i+1)/2, so that Eq. (37) is replaced by

Different from the original method developed by Brandt, the vector potential $\mathbf{A}$ was not involved during the derivation of above equations. Instead, $\mathbf{B}$ and $\mathbf{E}$ were used to avoid the unnecessary degree of freedom brought in by gauge selection. As a result, both arbitrary applied magnetic flux density $B_{a}(x, t)$ and applied current $I(t)$ are explicitly included in Eq. (48), making the method suitable for practical problems with strip geometry. As an example, the flux pump problem where complex applied field and applied transport current are preset conditions can be efficiently solved with this method [26].