††thanks: X.X. and Y.C. contributed equally to this work††thanks: X.X. and
Y.C. contributed equally to this work

# Scaling Law Breaking in Unequal-size Droplet Coalescence

Xi Xia  School of Mechanical Engineering, Shanghai Jiao Tong University,
Shanghai 200240, P. R. China  Yicheng Chi  School of Automotive and
Transportation Engineering, Shenzhen Polytechnic University, Shenzhen,
Guangdong 518055, P. R. China  Department of Mechanical Engineering, City
University of Hong Kong, Kowloon Tong 999077, Hong Kong  Peng Zhang
<EMAIL_ADDRESS>Department of Mechanical Engineering, City University of
Hong Kong, Kowloon Tong 999077, Hong Kong

###### Abstract

This Letter examines the coalescence of two unequal-size spherical liquid
droplets in the inviscid regime. We find that the liquid bridge evolution
exhibits a breaking from the classical 1/2 power-law scaling [Phys. Rev. Lett.
95, 164503 (2005)]. Employing an energy balance analysis, we attain a
theoretical model to collapse bridge evolution data of different droplet size
ratios. This model reveals an exponential dependence of the bridge’s radial
growth on time, which is intrinsically scaling-free owing to the asymmetric
movement of the liquid bridge.

In natural and industrial processes, the contact or impact of droplets Eggers
et al. (1999); Aarts et al. (2005); Yarin (2006); Thoraval et al. (2012); Tran
et al. (2013) could lead to the coalescence of liquid-gas interfaces, which is
crucial to the outcome or performance of the relevant applications. Extensive
research has been carried out to understand the most basic situation, the
momentumless coalescence of a pair of liquid droplets Thoroddsen et al.
(2005); Paulsen et al. (2011); Zhang and Law (2011); Kavehpour (2015). Early
studies Eggers et al. (1999); Duchemin et al. (2003) found that the radial
growth of the liquid bridge, which forms between the merging droplets,
satisfies certain scaling relations between the radius $R$ of the bridge and
time $t$. Later experimental Wu et al. (2004); Aarts et al. (2005); Thoroddsen
et al. (2005); Fezzaa and Wang (2008); Case (2009) and numerical Burton and
Taborek (2007); Pothier and Lewis (2012); Gross et al. (2013); Sprittles and
Shikhmurzaev (2012) works were able to confirm the existence of a $1/2$ power-
law scaling, i.e. $R\sim\sqrt{t}$, when the coalescence is in the inviscid (or
inertial) regime. On the other hand, the bridge evolution in the viscous
regime is better modeled by a linear scaling Aarts et al. (2005); Thoroddsen
et al. (2005); Yao et al. (2005); Burton and Taborek (2007); Paulsen et al.
(2011). The crossover (or transition) Burton and Taborek (2007) between the
viscous and inertial regimes has also attracted considerable interest, from
the discovery of a master curve for both regimes Paulsen et al. (2011);
Paulsen (2013) to the development of theoretical models justifying the
underlying universality Xia et al. (2019); Hack et al. (2020).

Previous research on binary droplet coalescence mainly revolves around two
equal-size droplets. However, less attention has been given to droplet
coalescence with size disparities, despite its higher relevance to reality.
Among the existing works involving the coalescence of unequal-size droplets
Anilkumar et al. (1991); Blanchette (2010); Liu et al. (2013); Sun et al.
(2015); Tang et al. (2016); Xia et al. (2017), the main focus was on the
effect of internal mixing facilitated by the break of symmetry. Regarding the
evolution of the liquid bridge, it is evident from our previous work Xia et
al. (2017) that the bridge surface shows an asymmetric growth–the bridge
interface becomes tilted as it expands out. Yet, little quantitative study
exits on the liquid bridge evolution of unequal-size droplet coalescence.

In this Letter, we study the liquid bridge evolution for two unequal-size
droplets merging in the inertial regime. We first conduct experiment to
resolve the bridge’s temporal variations for droplet pairs of various liquid
properties and size ratios. Then, we develop a theory to model the radial
growth of the liquid bridge, resulting in an exponential dependence of $R$ on
$t$ that breaks the classical $1/2$ power-law scaling relation.

Figure 1: (a) high-speed image of the sessile and pendant droplets and (b)
zoomed-in detail of the liquid bridge.

_Experimental approach $\&$ observations._—The present experiment employs the
classical sessile-pendant approach for droplet coalescence, similar to those
reported previously Thoroddsen et al. (2005); Aarts et al. (2005); Fezzaa and
Wang (2008); Case (2009); see the Supplemental Materials See Supplemental
Materials for a parameter list of all experimental cases.  for a schematic of
the setup. During each experimental run, a sessile droplet with a diameter of
1-5 mm is first generated by a syringe pump (Longer Precision Pump) and placed
on a super-hydrophobic surface, yielding a contact angle of $\sim 140^{\circ}$
and a near-spherical upper part as depicted in Fig. 1(a). Then, the syringe
pump generates a smaller-size pendant droplet (1-2 mm in diameter), which is
attached to the needle tip. Subsequently, the merging process is initiated by
actuating an automatic lifting platform (Winner Optics) which holds the super-
hydrophobic surface, slowly bringing the sessile droplet into contact with the
pendant droplet. The platform rises at a speed about 10 $\mu$m/s, which is
sufficiently small to ignore the gas-film flow disturbing the initial droplet
coalescence Case and Nagel (2008); Case (2009); Paulsen et al. (2011).

A high-speed camera (Photron SA-Z) integrated with a long-distance microscope
(Questar QM100) is used to capture time-resolved shadowgraph images of the
merging droplets. The camera operates at 150,000 frame per second (fps) with a
spatial resolution of 384$\times$256 pixels and a field of view of
2.04$\times$1.36 mm2, corresponding to a resolution of 5.3 $\mu$m/pixel. The
diameters of the small and large droplets, $D_{S}$ and $D_{L}$, are determined
at the initial moment prior to the coalescence, based on fitting an arc to
three arbitrarily-selected points on each droplet contour. The uncertainty
associated with the diameter measurement is estimated to be within $\pm 3\%$.
As such, the large-to-small droplet size ratio is defined as
$\Delta=D_{L}/D_{S}$, which varies from 1.0 to 5.0 in the present experiment.
Fig. 1(b) illustrates the zoomed-in detail of the liquid bridge, exhibiting a
distinct asymmetric bridge interface. So two characteristic radii of the
droplet bridge, $R_{S}$ and $R_{L}$, can be respectively defined as the radial
distances from the two points, $P_{1}$ and $P_{2}$, where the bridge interface
intersects the contours of the initial droplets (denoted by the white-dashed
lines), to the axis of symmetry. Then, the characteristic radius of the
circular bridge, $R$, is defined as $R=(R_{S}+R_{L})/2$. Note this definition
differs from the conventional one based on the minimum radial distance of the
bridge interface to the center axis, so it yields a slightly larger $R$.

To study the effect of varying liquid properties, i.e., density $\rho_{l}$,
dynamic viscosity $\mu_{l}$, and surface tension $\sigma$, we adopt water and
two aqueous glycerol solutions with 40 wt$\%$ and 60 wt$\%$ glycerol. The
three different liquids correspond to $\rho=$ 1000, 1100, and 1150
kg$\cdot$m-3, $\mu=$ 1.002, 3.630, and 10.80 mPa$\cdot$s, and $\sigma=$ 72.8,
70.0, and 66.0 mN$\cdot$m-1, respectively. The Ohnesorge number, defined as
$Oh=\mu(\rho\sigma D_{L})^{-1/2}$, varies within $10^{-3}$-$10^{-2}$. All
cases are characterized in terms of $\Delta$ and $Oh$, as summarized in the
Supplemental Materials See Supplemental Materials for a parameter list of all
experimental cases. .

The image sequences for representative droplet coalescence cases are presented
in Fig. 2(a-c), based on which the droplet interface contours corresponding to
the different snapshots are extracted and overlapped in Fig. 2(b-f). It is
seen that the liquid bridge displays a notable asymmetry for $\Delta>1$,
rendering an inclined bridge interface with $R_{L}>R_{S}$, which becomes
increasingly evident as the bridge evolves with time. As $Oh$ increases from
0.0016 to 0.0051, the bridge profiles in Fig. 2(f) follow the original droplet
contour more closely than those in Fig. 2(e), as the secondary deformation on
the bridge’s upper and lower surfaces tends to be inhibited. This can be
interpreted that the primary bridge movement gives rise to the development of
a capillary wave along the droplet surface, which can be damped by enhanced
viscosity. The experimental images for other cases can be found in the
Supplemental Materials See Supplemental Materials for a parameter list of all
experimental cases. .

Figure 2: Image sequences of the coalescence process of (a) equal-size water
droplets and (b-c) unequal-size droplets of water and 40 wt$\%$ aqueous
glycerol. (d-f) show the evolutions of the extracted interface contours
corresponding to (a-c).

_Geometric correlations._ —Inspired by the liquid bridge configuration
observed from the experiment, we present a geometric model for the bridge
surface between two unequal-size droplets in Fig. 3, based on the core
geometric variables of $R_{S}(t)$, $R_{L}(t)$, $\theta_{S}(t)$,
$\theta_{L}(t)$, and $S(t)$. In obtaining the geometric correlations, we
consider the following conditions.
(i) ‘Similar-size droplets’, meaning the two droplets are of similar sizes
such that $\Delta\sim O(1)$, $R_{L}/R_{S}\sim O(1)$, and
$\theta_{L}/\theta_{S}\sim O(1)$.
(ii) ‘Small bridge’, meaning that the characteristic radii of the bridge,
$R_{S}$ and $R_{L}$, are much smaller than the droplet diameters, i.e.,
$R_{S}/D_{S}\sim o(1)$ and $R_{L}/D_{L}\sim o(1)$. This condition is readily
satisfied during the early-stage coalescence.
(iii) ‘Equally-dividable bridge interface’, meaning there exists a principle
normal direction $\mathbf{n}_{p}$ dividing the bridge interface into similar-
shaped upper and lower sections, and the surface stress (including both
pressure difference $\Delta p$ and viscous stress $\mathbf{\tau}$) is also
equally distributed over the two sections. Given the present experimental
results, this condition can be considered a first approximation, although the
two sections divided by $\mathbf{n}_{p}$ might not be perfectly equal.

Based on Fig. 3, we can obtain the following geometric correlations,

$\begin{split}&\frac{R_{S}}{D_{S}}=\frac{\sin{(2\theta_{S})}}{2}=\theta_{S}-O(\theta_{S}^{3}),\\\
&\frac{R_{L}}{D_{L}}=\frac{\sin{(2\theta_{L})}}{2}=\theta_{L}-O(\theta_{L}^{3}).\end{split}$
(1)

Applying condition (ii), it is deduced that both $\theta_{S}$ and $\theta_{L}$
are of $o(1)$, $\theta_{S}\approx R_{S}/D_{S}$, and $\theta_{L}\approx
R_{L}/D_{L}$. It follows that the width of the bridge, $S$, defined as the
axial distance between $P_{1}$ and $P_{2}$, satisfies

$\begin{split}S&=R_{S}\tan{\theta_{S}}+R_{L}\tan{\theta_{L}}\\\
&=R_{S}\theta_{S}+R_{L}\theta_{L}+R_{S}O(\theta_{S}^{3})+R_{L}O(\theta_{L}^{3}).\end{split}$
(2)

Combining Eqs. LABEL:eq:m1 and 2, we can easily show that $S/D_{S}$ and
$S/D_{L}$ are of $O(\theta_{S}^{2})$ or $O(\theta_{L}^{2})$.

Figure 3: Schematic of the liquid bridge between two unequal-size merging
droplets.

Note that, compared with the equal-size situation, the unequal-size
coalescence is characterized by the tilted bridge interface and the misaligned
bridge movement from the radial direction, as illustrated in Fig. 3. The cause
of the tilted interface may be understood in a force balance analysis. Under
condition (iii), the principle normal direction $\mathbf{n}_{p}$ of the bridge
interface should align with the integral of the total surface stress across
the interface, which is balanced by the surface tension exerted at both ends
of the bridge interface. Given the two surface tension forces are identical to
each other, $\mathbf{n}_{p}$ must be in line with the angular bisector of the
two surface tension forces. Therefore, the tilting angle of the bridge
interface, $\theta$, can be simply expressed as

$\theta=\theta_{S}-\theta_{L}.$ (3)

Since both $\theta_{S}$ and $\theta_{L}\sim o(1)$, Eq. 3 implies that $\theta$
is also a $o(1)$ quantity. Thus, $\theta$ is treated as the ‘small parameter’
hereinafter. From Fig. 3, we further note another geometric relationship,
$R_{L}-R_{S}=S\tan{\theta}\approx S\theta$, which can be combined with Eqs.
LABEL:eq:m1 and 3 to yield

$\begin{split}\theta&\approx\left(\frac{1}{D_{S}}-\frac{1}{D_{L}}\right)R-\left(\frac{1}{D_{S}}+\frac{1}{D_{L}}\right)\frac{S\theta}{2}\\\
&\approx\left(\frac{1}{D_{S}}-\frac{1}{D_{L}}\right)R.\end{split}$ (4)

The second term on the right-hand side of Eq. 4 is negligible compared with
$\theta$ because $S/D_{S}\sim O(\theta_{S}^{2})$ and $S/D_{L}\sim
O(\theta_{L}^{2})$, as indicated from Eq. 2. Similarly, we can apply
$R_{L}-R_{S}\approx S\theta$ to Eqs. LABEL:eq:m1 and 2 to obtain

$\begin{split}S&\approx\left(\frac{1}{D_{S}}+\frac{1}{D_{L}}\right)R^{2}+\left(\frac{R_{L}+R}{D_{L}}-\frac{R_{S}+R}{D_{S}}\right)\frac{S\theta}{2}\\\
&\approx\left(\frac{1}{D_{S}}+\frac{1}{D_{L}}\right)R^{2}.\end{split}$ (5)

Again, the second term on the right-hand side of Eq. 5 is dropped out because
it is a $O(\theta^{2})S$ term.

On the other hand, given the equally-distributed viscous stress,
$\mathbf{\tau}=2\mu\mathbf{S}$ with $\mathbf{S}$ being the strain-rate tensor,
the liquid-side velocity must also be equally distributed over the bridge
interface, so the overall bridge movement also points in the direction of
$\mathbf{n}_{p}$. The physical implication here is significant. Consider the
very early stage of coalescence when $R_{S}\approx R_{L}$, Eq. 3 dictates that
$\theta$ is a positive value, which explains why the bridge movement is
inclined towards the smaller droplet from the beginning. It follows that the
velocity $V$ at which the bridge interface expands out is given by the
kinematic relationship, $\mathrm{d}R/\mathrm{d}t=V\cos{\theta}$. After Taylor
expansion, it takes the form

$\frac{\mathrm{d}R}{\mathrm{d}t}\approx V\left(1-\frac{\theta^{2}}{2}\right).$
(6)

_Energy balance analysis and solution._ —Now, considering the energy balance
during the coalescence process, the movement of the liquid entrained by the
bridge is driven by the rapid discharge of the surface energy, which can be
expressed as

$\Delta E_{s}+\Delta E_{k}\approx 0,$ (7)

where $\Delta E_{s}$ and $\Delta E_{k}$ represent the changes in surface and
kinetic energies, respectively, from the initial state before coalescence.
Note Eq. 7 requires the viscous dissipation being negligible compared to the
change in inertia energy. Based on the geometries in Fig. 3, we have

$\begin{split}\Delta E_{s}&=-\pi
D_{S}^{2}(\theta_{S}^{2}+O(\theta_{S}^{4}))\sigma-\pi
D_{L}^{2}(\theta_{L}^{2}+O(\theta_{L}^{4}))\sigma\\\ &+\pi
S(R_{S}+R_{L})[1+O((\theta_{S}+\theta_{L})^{2})]\sigma,\end{split}$ (8)

where $\sigma$ is the surface tension. Since $\theta_{S}$ and $\theta_{L}\sim
o(1)$, we can combine Eqs. LABEL:eq:m1 and 2 to obtain $S\approx
R_{S}\theta_{S}+R_{L}\theta_{L}\approx
D_{S}\theta_{S}^{2}+D_{L}\theta_{L}^{2}$. Thus, the third term on the right-
hand side of Eq. 8 has the leading order of $O(\theta_{S}^{3})D_{S}^{2}\sigma$
or $O(\theta_{L}^{3})D_{L}^{2}\sigma$, which is one order smaller than the
first two terms and can be dropped out. Next, $\Delta E_{k}$ can be estimated
as

$\Delta E_{k}\approx\frac{C}{2}\pi\rho_{l}R^{2}SV^{2},$ (9)

where $\rho_{l}$ is the liquid density and $C$ is a prefactor related to the
volume and velocity distribution of the moving fluid entrained with the
bridge.

With $S$ given by Eq. 5, we can plug Eqs. 8 and 9 into Eq. 7 and balance the
leading order terms to obtain

$\frac{C}{2}\rho_{l}R^{4}\left(\frac{1}{D_{S}}+\frac{1}{D_{L}}\right)V^{2}\approx
2R^{2}\sigma.$ (10)

Further substituting $V$ using Eqs. 6 and 4, we have

$\frac{1}{1-\beta
R^{2}/2}\frac{\mathrm{d}R}{\mathrm{d}t}\approx\frac{\gamma}{R}.$ (11)

In simplifying Eq. 11, we define $\beta$ and $\gamma$ as

$\begin{split}\beta=\left(\frac{1}{D_{S}}-\frac{1}{D_{L}}\right)^{2}\text{ and
}\gamma=\left[\frac{4\sigma
D_{S}}{C\rho_{l}(1+1/\Delta)}\right]^{1/2}.\end{split}$ (12)

Eq. 11 has the solution:

$R^{2}\sim\frac{1-e^{-\beta\gamma t}}{\beta},$ (13)

where the initial condition, $R(t=0)=0$, has been used in the derivation. Eq.
13 can further take the nondimensional form of

$R^{*2}\sim 1-e^{-t^{*}},$ (14)

where $R^{*}=R\beta^{1/2}$ and $t^{*}=t\beta\gamma$. Interestingly, the bridge
evolution governed by Eq. 13 or 14 no longer has a power-law scaling between
$R$ and $t$. However, by letting $\Delta\rightarrow 1$, we have
$\beta\rightarrow 0$ and $R^{2}\sim[1-(1-\beta\gamma t)]/\beta=\gamma t$. This
means that the present model is able to recover the inviscid scaling law of
$R\sim t^{1/2}$ in the equal-size limit.

_Discussion._ —It is worth discussing the origin of this scaling-free
solution. Mathematically, the exponential dependence of $R$ on $t$ arises from
the $\theta^{2}$ term in Eq. 6. If $\theta=0$ is set in Eq. 6, the above
analysis would give the exact scaling relation of $R\sim t^{1/2}$. Thus, we
can infer that the breaking of the power-law scaling results from the
misaligned movement of the liquid bridge from the radial direction.

Next, we check the scaling breaking behavior for the present experimental data
of various $\Delta$ and $Oh$. In Fig. 4, the classical scaling law of $R\sim
t^{1/2}$ is evaluated by plotting the data in the parameter space
$[(8\sigma/\rho D_{S}^{3})^{1/2}t,(2R/D_{S})^{2}]$ Duchemin et al. (2003);
Aarts et al. (2005). While the overall scaling of $R^{2}\sim t$ is still valid
for most cases, there is an apparent upward drift of large-$\Delta$ data. And
a larger $\Delta$ tends to cause a larger shift from the $\Delta=1.0$ line.
This can be understood that the presence of a larger droplet enhances the
expansion speed of the bridge interface, corresponding to a larger prefactor
of the power-law scaling. Furthermore, the cases with $\Delta>3$ display a
gradual deflection from $R^{2}\sim t$ as time proceeds and the liquid bridge
becomes more asymmetric, indicating a diversion from the power-law scaling.

Figure 4: Effect of $\Delta$ on the scaling of $R^{2}\sim t$. See Supplemental
Materials See Supplemental Materials for a parameter list of all experimental
cases.  for detailed experimental parameters. Figure 5: Validation of Eq. 14
against experimental data.

To evaluate our theory, Fig. 5 shows the validation of Eq. 14 based on the
same set of data in Fig. 4. We can confirm the collapse of non-unity-$\Delta$
data onto a single line given by Eq. 14 with $C=1.5$ ($C$ is obtained by
fitting Eq. 14 with a unity prefactor). Here, the effect of $\Delta$ is
assimilated into the time and length scales. In particular, a larger $\Delta$
generally corresponds to higher $R^{*}$, as $R$ is scaled by $\beta^{-1/2}$
that is highly sensitive to $\Delta$. This result also demonstrates that the
present model has a negligible dependence on viscosity for $Oh$ varying
between $10^{-3}$ and $10^{-2}$. The theoretical line in Fig. 5 remains almost
linear when $R^{*2}<10^{-1}$, while it exhibits notable deflection or
deviation from the power-law scaling when $R^{*2}>10^{-1}$. Given the
correlation between $\Delta$ and $R^{*}$, this means that the scaling-law
breaking is manifest when $\Delta$ is greater than 3 or so, which explains the
observed scaling diversion in Fig. 4.

In summary, this Letter reports the breaking of the $1/2$ power-law scaling
for the liquid bridge evolution during the coalescence of unequal-size
droplets. By accounting for the asymmetric motion of the liquid bridge
interface, we attain a new solution for the bridge evolution in the inviscid
regime based on balancing the changes in surface and kinetic energies. The
derived model shows distinct scaling-free feature and well resolves the
droplet coalescence of non-unity size ratios.

We acknowledge financial support from the National Natural Science Foundation
of China (No. 52176134 and 12072194), the Research Grants Council of the Hong
Kong Special Administrative Region, China (Project No. CityU 15218820), and
the APRC-CityU New Research Initiatives/Infrastructure Support from Central of
City University of Hong Kong (No. 9610601). We thank Dr. Tao Yang and Jiayue
Han for their assistance in data processing.

## References

  * Eggers et al. (1999) J. Eggers, J. R. Lister, and H. A. Stone, J. Fluid Mech. 401, 293 (1999).
  * Aarts et al. (2005) D. G. A. L. Aarts, H. N. W. Lekkerkerker, H. Guo, G. H. Wegdam, and D. Bonn, Phys. Rev. Lett. 95, 164503 (2005).
  * Yarin (2006) A. L. Yarin, Ann. Rev. Fluid Mech. 38, 159 (2006).
  * Thoraval et al. (2012) M.-J. Thoraval, K. Takehara, T. G. Etoh, S. Popinet, P. Ray, C. Josserand, S. Zaleski, and S. T. Thoroddsen, Phys. Rev. Lett. 108, 264506 (2012).
  * Tran et al. (2013) T. Tran, H. de Maleprade, C. Sun, and D. Lohse, J. Fluid Mech. 726, R3 (2013).
  * Thoroddsen et al. (2005) S. T. Thoroddsen, K. Takehara, and T. G. Etoh, J. Fluid Mech. 527, 85 (2005).
  * Paulsen et al. (2011) J. D. Paulsen, J. C. Burton, and S. R. Nagel, Phys. Rev. Lett. 106, 114501 (2011).
  * Zhang and Law (2011) P. Zhang and C. K. Law, Phys. Fluids 23, 042102 (2011).
  * Kavehpour (2015) H. P. Kavehpour, Ann. Rev. Fluid Mech. 47, 245 (2015).
  * Duchemin et al. (2003) L. Duchemin, J. Eggers, and C. Josseran, J. Fluid Mech. 487, 167 (2003).
  * Wu et al. (2004) M. Wu, T. Cubaud, and C.-M. Ho, Phys. Fluids 16, L51 (2004).
  * Fezzaa and Wang (2008) K. Fezzaa and Y. Wang, Phys. Rev. Lett. 100, 104501 (2008).
  * Case (2009) S. C. Case, Phys. Rev. E 79, 026307 (2009).
  * Burton and Taborek (2007) J. C. Burton and P. Taborek, Phys. Rev. Lett. 98, 224502 (2007).
  * Pothier and Lewis (2012) J. C. Pothier and L. J. Lewis, Phys. Rev. B 85, 115447 (2012).
  * Gross et al. (2013) M. Gross, I. Steinbach, D. Raabe, and F. Varnik, Phys. Fluids 25, 052101 (2013).
  * Sprittles and Shikhmurzaev (2012) J. E. Sprittles and Y. D. Shikhmurzaev, Phys. Fluids 24, 122105 (2012).
  * Yao et al. (2005) W. Yao, H. J. Maris, P. Pennington, and G. M. Seidel, Phys. Rev. E 71, 016309 (2005).
  * Paulsen (2013) J. D. Paulsen, Phys. Rev. E 88, 063010 (2013).
  * Xia et al. (2019) X. Xia, C. He, and P. Zhang, Proc. Natl. Acad. Sci. U.S.A. 116, 23467 (2019).
  * Hack et al. (2020) M. A. Hack, W. Tewes, Q. Xie, C. Datt, K. Harth, J. Harting, and J. H. Snoeijer, Phys. Rev. Lett. 124, 194502 (2020).
  * Anilkumar et al. (1991) A. V. Anilkumar, C. P. Lee, and T. G. Wang, Phys. Fluids A 3, 2587 (1991).
  * Blanchette (2010) F. Blanchette, Phys. Rev. Lett. 105, 074501 (2010).
  * Liu et al. (2013) D. Liu, P. Zhang, C. K. Law, and Y. Guo, Int. J. Heat Mass Transfer 57, 421 (2013).
  * Sun et al. (2015) K. Sun, T. Wang, P. Zhang, and C. K. Law, Phys. Rev. E 91, 023009 (2015).
  * Tang et al. (2016) C. Tang, J. Zhao, P. Zhang, C. K. Law, and Z. Huang, J. Fluid Mech. 795, 671 (2016).
  * Xia et al. (2017) X. Xia, C. He, D. Yu, J. Zhao, and P. Zhang, Phys. Rev. Fluids 2, 113607 (2017).
  * Case and Nagel (2008) S. C. Case and S. R. Nagel, Phys. Rev. Lett. 100, 084503 (2008).
  * (29) See Supplemental Materials for a parameter list of all experimental cases.


#  V889 Her: abrupt changes in the magnetic field or differential rotation?

Teemu Willamo  Thomas Hackman  Department of Physics, P.O. Box 64, FI-00014
University of Helsinki, Finland  Jyri J. Lehtinen  Finnish Centre for
Astronomy with ESO (FINCA), University of Turku, Vesilinnantie 5, FI-20014
University of Turku, Finland  Department of Physics, P.O. Box 64, FI-00014
University of Helsinki, Finland  Maarit J. Korpi-Lagg  Department of Computer
Science, Aalto University, PO Box 15400, FI-00076 Aalto, Finland  Max Planck
Institute for Solar System Research, Justus-von-Liebig-Weg 3, D-37077
Göttingen, Germany  Oleg Kochukhov  Department of Physics and Astronomy,
Uppsala University, Box 516, 751 20 Uppsala, Sweden

###### Abstract

We have applied Zeeman-Doppler imaging (ZDI) to an extensive
spectropolarimetric HARPSpol data set of the magnetically active young solar
analogue V889 Her, covering 35 spectra obtained during six nights in May 2011.
The data set allows us to study Stokes V profiles of the star at almost
identical rotational phases, separated by one or more stellar rotations. We
use these data to study if the line profiles evolve from one rotation to the
next, and find that some evolution does indeed occur. We consider two possible
explanations for this: abrupt changes in the large-scale magnetic field or
differential rotation. We find it quite difficult to distinguish between the
two alternatives using ZDI alone. A strong differential rotation could,
however, explain the changes in the line profiles, so we conclude that it must
be present, and the abrupt magnetic field evolution is left uncertain.
Commonly, rapidly rotating stars are assumed to have only weak differential
rotation. If the strong differential rotation of V889 Her is indeed present,
as has been found in other studies as well, it could indicate that the
theoretical and numerical results of differential rotation still need to be
revised. The rapid changes that may have occurred in the magnetic field
indicate that one should be quite cautious when interpreting ZDI maps
constructed from data over long time intervals.

###### keywords:

stars: activity — stars: individual (V889 Her) — stars: magnetic field —
stars: solar-type

##  1 Introduction

Zeeman-Doppler imaging (ZDI) is a powerful tool used to indirectly map the
magnetic field at the stellar surface (e.g. Brown et al., 1991; Kochukhov,
2016). This is a further development of Doppler imaging (DI; Vogt et al.,
1987; Piskunov et al., 1990), which is used to map the surface temperature
distribution. The magnetic field is mapped by inverting a series of Stokes I
and V profiles (and in some cases also Q and U profiles) of absorption lines,
distributed over the stellar rotation phases, thus finding the magnetic field
and brightness distribution which reproduces the observed rotational changes
in the line profiles. To map the whole visible stellar surface, a good
coverage over the rotational phases is required. Depending on the stellar
rotation period, this requires observations over many nights, especally for
more slowly rotating stars or stars with a rotation period close to a small
multiple of one day. In order to produce a sensible ZDI map, one has to assume
that the magnetic field does not evolve during these observations. Similarly,
with DI, it is assumed that the spot configuration stays constant through the
observations. However, Strassmeier and Bartus (2000) showed that there has
been some rapid spot evolution happening in the RS CVn type star V711 Tau.
Strassmeier et al. (2003a) analyzed two DI maps from two consecutive rotations
of the K giant V1192 Ori, and noted some differences, which they concluded to
be due to anti-solar differential rotation. Donati et al. (2016, 2017) noted
differences in the magnetic topology of two data sets of the T Tauri type star
V830 Tau, which were separated by approximately one month, and were explained
with differential rotation. Also Lehtinen et al. (2022) noted some evolution
of the Stokes V profiles of the young solar analogue LQ Hya in two data sets
with the duration of 9 and 11 consecutive nights.

The rate of changes happening in the magnetic field corresponds to the size of
the magnetic structures. Large structures generally evolve slowly, while small
structures can evolve more rapidly (Giles et al., 2017). Small-scale fields,
however, are difficult to detect with ZDI, since opposite polarities of nearby
regions cancel each other out. At worst, ZDI is able to recover less than 1
$\%$ of the total magnetic field energy (Kochukhov et al., 2020).

The object of this study, V889 Her, is a magnetically active ($\log
R^{\prime}_{\mathrm{HK}}=-4.175$; Lehtinen et al., 2016), young solar
analogue, of spectral class G2 V (Montes et al., 2001) and an estimated age of
30-50 million years (Strassmeier et al., 2003b). There are a multitude of
published DI and ZDI maps of V889 Her, most recently by Willamo et al. (2019).
They published 19 temperature maps of V889 Her, where the main feature was a
dark polar spot present in all of them. Previous DI maps have been published
by Strassmeier et al. (2003b), Järvinen et al. (2008), Huber et al. (2009) and
Frasca et al. (2010), and ZDI maps by Marsden et al. (2006) and Jeffers and
Donati (2008), who both found also a strong differential rotation to be
present. Huber et al. (2009) found no evidence for differential rotation, and
Kővári et al. (2011) found a differential rotation that was an order of
magnitude weaker than that found by Marsden et al. (2006) and Jeffers and
Donati (2008).

According to both theoretical (Kitchatinov and Rüdiger, 1999) and numerical
(e.g. Viviani et al., 2018) studies, rapidly rotating stars should not have
very strong differential rotation relative to their rotation period. Many
observational studies also agree with this (e.g. Reiners and Schmitt, 2003;
Lehtinen et al., 2016), although, as seen with the conflicting results of V889
Her, the observational results are not as clear. In the study of Reinhold and
Gizon (2015) we can see that the differential rotation relative to the
rotation rate decreases for rapid rotators, although the absolute rotational
shear actually slightly increases. All in all, differential rotation has
proven to be a difficult parameter to measure accurately from real stars (see
also e.g. Hackman et al., 2019). Since differential rotation is a large
unknown factor, it is unclear to what extent it can explain apparent evolution
of spectral line profiles, and when these apparent changes correspond to true
evolution of the magnetic field or spot structures.

##  2 Data

We analyze data of V889 Her obtained with the HARPS-pol high-resolution
spectropolarimeter (Piskunov et al., 2011) at the ESO 3.6m telescope at La
Silla, Chile, retrieved from the ESO archives. The observations are from May
15–20, 2011. They were reduced with the REDUCE package (Piskunov and Valenti,
2002). Spread over a period of 6 nights and containing a total of 35
individual spectra, this is a data set of remarkably high quality, in terms of
the signal-to-noise ratio and phase coverage. See Table 1 for more details.
The phase coverage was enhanced by combining only two sub-exposures, instead
of the standard four, for the production of LSD profiles (see next Section).
Due to this, no null profiles are available for the spectra.

To calculate the rotation phase $\phi$ of the star we have used the ephemeris:

$\mathrm{HJD}_{\phi=0}=2449950.550+P_{\mathrm{rot}}\times E.$ (1)

The rotation period $P_{\mathrm{rot}}=1.3300\pm 0.0031$ d was derived using
photometry from Willamo et al. (2019), including only the 28 data points
between Julian Dates 2455685 and 2455710, corresponding to May 2011. This way
we get the best period for this particular data set, and are not affected by
fluctuations in the period due to differential rotation or other factors. We
determined the best period with the Three Stage Period Analysis method (Jetsu
and Pelt, 1999).

Table 1: Observation log for V889 Her. #Night numbers the observing nights as they are referred to in the text. $S/N$ is the signal-to-noise ratio of the Stokes LSD V-profile. #Night | Date (UT) | HJD-2455000 | $\phi$ | $S/N$
---|---|---|---|---
1 | 2011 05 15 5:49 | 696.743 | 0.446 | 41833
1 | 2011 05 15 6:17 | 696.763 | 0.461 | 40821
1 | 2011 05 15 7:51 | 696.828 | 0.510 | 43419
1 | 2011 05 15 8:19 | 696.847 | 0.524 | 44942
1 | 2011 05 15 9:42 | 696.905 | 0.568 | 37110
1 | 2011 05 15 10:04 | 696.920 | 0.579 | 36540
1 | 2011 05 15 10:27 | 696.936 | 0.591 | 36018
2 | 2011 05 16 5:50 | 697.744 | 0.198 | 42084
2 | 2011 05 16 6:18 | 697.763 | 0.213 | 42375
2 | 2011 05 16 7:50 | 697.827 | 0.261 | 42257
2 | 2011 05 16 8:18 | 697.847 | 0.276 | 44445
2 | 2011 05 16 9:47 | 697.908 | 0.322 | 42123
2 | 2011 05 16 10:15 | 697.928 | 0.337 | 42839
3 | 2011 05 17 5:19 | 698.722 | 0.934 | 41581
3 | 2011 05 17 5:47 | 698.742 | 0.949 | 44364
3 | 2011 05 17 8:18 | 698.847 | 0.028 | 39113
3 | 2011 05 17 8:40 | 698.862 | 0.039 | 37740
3 | 2011 05 17 10:05 | 698.921 | 0.083 | 28415
3 | 2011 05 17 10:26 | 698.936 | 0.095 | 27794
4 | 2011 05 18 5:23 | 699.725 | 0.688 | 33248
4 | 2011 05 18 5:51 | 699.745 | 0.703 | 37410
4 | 2011 05 18 8:17 | 699.846 | 0.779 | 37978
4 | 2011 05 18 8:38 | 699.861 | 0.790 | 38071
4 | 2011 05 18 9:00 | 699.876 | 0.802 | 36297
5 | 2011 05 19 4:54 | 700.705 | 0.425 | 34808
5 | 2011 05 19 5:21 | 700.724 | 0.439 | 31547
5 | 2011 05 19 7:55 | 700.831 | 0.520 | 35895
5 | 2011 05 19 8:23 | 700.850 | 0.534 | 32362
5 | 2011 05 19 8:52 | 700.870 | 0.549 | 34577
6 | 2011 05 20 5:19 | 701.722 | 0.189 | 24631
6 | 2011 05 20 5:47 | 701.742 | 0.205 | 26713
6 | 2011 05 20 7:17 | 701.804 | 0.251 | 23753
6 | 2011 05 20 7:45 | 701.824 | 0.266 | 18648
6 | 2011 05 20 9:31 | 701.897 | 0.321 | 21406
6 | 2011 05 20 9:59 | 701.917 | 0.336 | 25625

##  3 Zeeman-Doppler imaging

As the polarization signal in individual spectral lines, caused by the
magnetic field, is too weak to be reliably detected, we use the Least-Square
Deconvolution (LSD) method to enhance the signal. The method works by
combining thousands of spectral lines into a mean LSD profile (Kochukhov et
al., 2010). The line parameters are extracted from the VALD
database111http://vald.astro.uu.se/ (Piskunov et al., 1995; Kupka et al.,
1999). In the final profiles, 6545 lines between 3900-7100 $\mathrm{\AA}$ were
used, with a mean wavelength 5111.589 $\mathrm{\AA}$ and a mean Landé
$g$-factor 1.215.

We use the code inversLSD, developed by Kochukhov et al. (2014), to invert the
ZDI maps. The stellar surface is divided into 1876 elements with approximately
equal sizes, and the radial, meridional and azimuthal magnetic field
components and the brightness is calculated for each. Instead of modeling
directly the magnetic field components, the magnetic field is represented with
spherical harmonic expansions , where higher order harmonics represent the
magnetic field on smaller scales. For the number of spherical harmonic
expansions, Hackman et al. (2016) found that increasing maximum angular degree
$\ell_{\mathrm{max}}$ to more than 5 did not significantly improve the fit to
the Stokes V profiles. Also Lehtinen et al. (2022) noted, that although they
used $\ell_{\mathrm{max}}=20$, most of the magnetic energy was found in
harmonics below $\ell\leq 6$. Both these studies analyzed stars with somewhat
different values for $v\sin i$ than V889 Her, but still of the same order of
magnitude, so these results should be fairly comparable to V889 Her. Similarly
to Willamo et al. (2022), we have limited the solution to
$\ell_{\mathrm{max}}=10$ harmonics.

The magnetic field and brightness inversions are performed simultaneously, so
that brightness inhomogeneities are taken into account in the magnetic field
inversion. This is necessary, since a lower brightness will lead to a weaker
polarization signal, which would be difficult to tell apart from a weaker
magnetic field, if the inversions were done independently of each other. For
the magnetic field, high-order modes are suppressed by a regularization
function, while for the brightness inversion, Tikhonov regularization is used
to constrain the solution to distributions with a minimal gradient. Smoother
variations are thus preferred over sharper ones. The regularization and other
parameters used in the inversion, as well as the stellar parameters, are shown
in Table 2. A more thorough description of the regularization can be found in
Rosén et al. (2016).

In our previous study of V889 Her (Willamo et al., 2019), the inclination $i$
was set to 60∘, following Marsden et al. (2006). However, from the light curve
recently obtained by the TESS space telescope222Available at
https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html, V889 Her seems
to have at least one transiting exoplanet, which would indicate that the
inclination is close to 90∘, if the rotation of the star is aligned with the
planetary disk. Thus we set the inclination to a higher value, $i=70^{\circ}$,
estimated with the equation:

$\sin i=\frac{P_{\mathrm{rot}}v\sin i}{2\pi R},$ (2)

where $P_{\mathrm{rot}}$ is the rotation period, $v\sin i$ the projected
rotation velocity, and $R$ the stellar radius. References for these parameters
are given in Table 2.

Table 2: Stellar parameters of V889 Her and parameters used in the ZDI inversions. Parameter | Value | Reference
---|---|---
$P_{\mathrm{rot}}$ | 1.3300 d | This study
$v\sin i$ | 38.5 km/s | Willamo et al. (2019)
$i$ | $70^{\circ}$ | This study
$R$ | 1.09 $R_{\odot}$ | Strassmeier et al. (2003b)
Microturbulence | 1.6 km/s | Järvinen et al. (2008)
Macroturbulence | 3.0 km/s | Strassmeier et al. (2003b)
Regularization (brightness) | $5\times 10^{-9}$ | …
Regularization (magnetic field) | $8\times 10^{-13}$ | …
Effective Landé $g$ factor | 1.215 | …
Mean $\lambda$ | 5630.0 Å | …

##  4 Differential rotation

The latitudinal differential rotation parameter $\alpha$ is defined as:

$\alpha=\frac{\Omega_{\mathrm{eq}}-\Omega_{\mathrm{pol}}}{\Omega_{\mathrm{eq}}},$
(3)

and the absolute latitudinal rotational shear is:

$\Delta\Omega=\Omega_{\mathrm{eq}}-\Omega_{\mathrm{pol}},$ (4)

where $\Omega_{\mathrm{eq}}$ and $\Omega_{\mathrm{pol}}$ are the angular
velocities at the equator and at the poles, respectively. Published values for
the differential rotation of V889 Her are $\Delta\Omega=0.402\pm 0.044$ rad/d
(Marsden et al., 2006), $\Delta\Omega=0.52\pm 0.04$ rad/d for Stokes I data
and $\Delta\Omega=0.47\pm 0.04$ rad/d for Stokes V data (Jeffers and Donati,
2008), and $\Delta\Omega=0.042$ rad/d (Kővári et al., 2011). Järvinen et al.
(2008) found a weaker differential rotation than Marsden et al. (2006),
although they did not derive a quantitative value for it, and Huber et al.
(2009) found no evidence for differential rotation. In Willamo et al. (2019),
V889 Her was assumed to have no differential rotation.

The differential rotation law is assumed to follow the solar differential
rotation law, with the higher order terms neglected. The dependence of the
surface rotation rate on the latitude $\theta$ will thus be:

$\Omega(\theta)=\Omega_{\mathrm{eq}}(1-\alpha\sin^{2}(\theta)).$ (5)

The differential rotation does not need to necessarily follow the solar
differential rotation law, but it is assumed, since we have no information to
make reliable guesses about alternative differential rotation laws.

With different values of $\alpha$, the rotation period $P_{\mathrm{rot}}$ must
also be changed. Since $P_{\mathrm{rot}}$ in the model is measured at the
equator, but the observational value comes from the latitude of the dominating
spot structure, the optimal value of $P_{\mathrm{rot}}$ is expected to
decrease with increasing $\alpha$. In practice, the best values for $\alpha$
and $P_{\mathrm{rot}}$ are found by trying a grid of different values, and
finding the minimum for the deviation between model and observations.

##  5 Stokes V profiles

The data set analyzed here contains 35 spectra. A typical (Z)DI data set
contains perhaps a little more than 10 spectra (e.g. Willamo et al., 2019;
Hackman et al., 2019; Lehtinen et al., 2022; Willamo et al., 2022), which, if
spread evenly over the rotational phases, already provide a sufficient phase
coverage. There are also studies with similar amounts of phases as here,
though, for example of V830 Tau (Donati et al., 2017). Our data set of V889
Her also contains some spectra at almost identical phases, which are separated
by one or a few rotations. These are essential for the aim of this study,
since, if no magnetic field evolution or differential rotation is present, the
surface structure seen to the observer, and thus also the line profiles, from
these phases should be identical. These multiple observations of the same
phase can thus be used as a rough check of the assumption of the stability of
the magnetic field over the time span of the observations. They are, however,
not very useful for improving the phase coverage, and are thus commonly
avoided in the observations due to time limitations, especially when one tries
to get sufficient observations of as many stars as possible.

Fig. 1 shows the Stokes LSD I and V profiles of the spectra, and the model
fitted to the data, both separately for the first three nights and last three
nights. There is especially one pair of profiles at $\phi=0.337$ and
$\phi=0.336$, where the Stokes V profiles are clearly not identical, although
separated only by $\Delta\phi=0.001$ (shown together in the lower panel of
Fig. 1). The difference between these is evidence that the star cannot have
rotated as a solid body with a constant magnetic field structure during these
observations, assuming that the rotation period we use is correct. Also phases
around $\phi=0.2$ (0.198 and 0.213 in the first three nights, 0.189 and 0.205
in the last three nights) and $\phi=0.5$ (0.510 and 0.524 in the first three
nights, 0.520 in the last three nights) are notably different in the figures.

Active stars, especially M dwarfs, are often showing flaring activity. Whereas
flares can affect Stokes I profiles, they should not significantly affect the
polarized Stokes V profiles (e.g. Donati et al., 2017). As a G-type star, V889
Her is not expected to show very high levels of flaring activity. A flare
would also raise the continuum level of the I profile, which would be seen as
a significant weakening of the line. This is not seen in any of the spectra.
Thus, flares cannot explain the discrepancies seen at simultaneous rotational
phases, as there are multiple examples of differences in the Stokes V profiles
of closeby phases.

Figure 1: Upper panel: the Stokes LSD profiles for the three first nights
(left) and the three last nights (right). The black lines are the observed LSD
profiles, and the red lines are fits made to these. The rotational phase is
indicated for each profile. Lower panel: V profiles for the two phases closest
to each other, where the difference is most notable. The phase $\phi=0.337$ is
part of the map for nights 1-3, while the phase $\phi=0.336$ is part of the
map for nights 4-6, in Fig. 2. This indicates that either a strong
differential rotation must be present, or there has been some rapid evolution
of the magnetic field.

##  6 Zeeman-Doppler maps

In this section we present the ZDI maps constructed from our data. We first
study the case where the star is assumed to rotate as a solid body with no
differential rotation and the magnetic field is allowed to evolve between two
subsets of the observations. After this we require the magnetic field to stay
constant through the observations, and allow for differential rotation. These
alternatives are then compared to each other, using the deviation between the
model and observations:

$\sigma=\frac{\sum_{\phi,\lambda}\omega_{\phi,\lambda}(I_{\mathrm{obs},\phi,\lambda}-I_{\mathrm{mod},\phi,\lambda})^{2}}{n_{\phi}n_{\lambda}},$
(6)

where $I_{\mathrm{obs},\phi,\lambda}$ and $I_{\mathrm{mod},\phi,\lambda}$ are
the observed and modeled intensities at each phase $\phi$ and wavelength
$\lambda$, $\omega_{\phi,\lambda}$ is the weight of each data point, which is
proportional to the square of the $S/N$ of the spectrum, and $n_{\phi}$ and
$n_{\lambda}$ are the number of phases and wavelength points. $\sigma$ is
shown for each model in Table 3, where we see that the static model with no
magnetic field evolution or differential rotation clearly has the highest
deviation.

###  6.1 Case 1: Abruptly evolving magnetic field

First we perform the analysis without differential rotation, but allowing the
magnetic field to evolve during the observations.

As it was clear from Fig. 1 that without differential rotation the magnetic
field had to evolve during the observations, we divided our observation set in
two, to get two subsets containing three consecutive nights each, labeled as
‘1-3’ and ‘4-6’. We also divided a third subset from the data, from which we
removed the spectra of the first and last night, thus containing four nights
from the middle of the observation set, labeled as ‘2-5’. It might be
questionable whether the magnetic field stays constant even for the period of
three or four nights, but by producing separate maps for these subsets, we
show that the resulting fields are indeed different.

If we assume one spectrum to cover 10 % of the rotation period, then all six
nights combined have a phase coverage $f_{\phi}=0.967$. With the division of
the observations to three subsets, we have phase coverages
$f_{\phi,1-3}=0.744$ for the first set, $f_{\phi,2-5}=0.925$ for the middle
set and $f_{\phi,4-6}=0.674$ for the last set, all being acceptable phase
coverages on their own. We used the same inversion parameters for each subset.

The three resulting ZDI maps are shown in Fig. 2, separately for the radial,
meridional and azimuthal field components, and the brightness. We see that the
subsets 1-3 and 2-5 are very similar, but differ from the subset 4-6, where
the magnetic field strength is reduced significantly. If this is the case,
then the magnetic field evolution would have happened quite abruptly during
the last nights.

The deviation $\sigma$ for each of these maps is shown in Table 3. To study
whether an abruptly changing magnetic field or differential rotation is a
better explanation to the data, these values should be compared to the one
resulting from the analysis involving differential rotation, which is studied
in the next section.

###  6.2 Case 2: Differential rotation

Now we test how well the data can be reconstructed when differential rotation
is allowed. This will stretch the magnetic field, which is assumed to be
frozen into the stellar surface, but it is otherwise kept constant. A map of
the deviation $\sigma$ for different values of $\alpha$ and $P_{\mathrm{rot}}$
is shown in Fig. 3. The minimum for $\sigma$ is found with parameters
$\alpha=0.075$333This corresponds to $\Delta\Omega=0.36$ rad/d, slightly less
than the values derived by Marsden et al. (2006) and Jeffers and Donati
(2008). and $P_{\mathrm{rot}}=1.31$ d, with $\sigma=2.5298\times 10^{-5}$. The
ZDI map corresponding to these best fit parameters is shown in Fig. 4. The
corresponding LSD profiles are shown in Fig. 5. When looking at the profiles,
the fits seem fairly good, which could mean that a strong differential
rotation can explain the observed apparent line profile evolution even without
significant abrupt changes in the magnetic field.

The uncertainty of the differential rotation is difficult to estimate
quantitatively. In Fig. 6 we show how $\sigma$ changes as a function of
$\alpha$, where the best $P_{\mathrm{rot}}$ value is always chosen for each
$\alpha$. The minimum of $\sigma$ is at $\alpha=0.075$, but with a larger
$\alpha$, the value of $\sigma$ increases only a little. We did thus not
estimate any formal error bars for $\alpha$. This demonstrates further the
difficulties in accurately estimating the amount of differential rotation,
even for very high-quality data.

Figure 2: Magnetic field and brightness maps for nights 1-3 (up), 2-5 (centre)
and 4-6 (down), without differential rotation. The vertical lines correspond
to phases used for each map. The horizontal line indicates the limit of the
visible surface due to the inclination. Figure 3: Map of the deviation
$\sigma$ for differential rotation parameters of V889 Her. The colour scaling
corresponds to the deviation between the model and observations in the
respective ZDI map, with the unit being $10^{-5}$. The best fit is indicated
with a cross.

From Table 3 we see that the case with differential rotation is a better fit
to the data than the subset 4-6 without differential rotation (note, however,
that night 6 has slightly lower S/N than the other nights, which could
somewhat increase the value of $\sigma$ for data sets where this night is
included), but subsets 1-3 and 2-5 fit the data better than the model with
differential rotation. A clear distinction between the cases is thus quite
difficult.

The results above were derived from only Stokes V data. We performed a similar
study on I data as well. The results from the I data, however, indicate a very
strong anti-solar-like differential rotation (negative $\alpha$). This is
probably artificially caused by the cool polar spot, which has a similar
influence on the mean line profile as anti-solar differential rotation would
have (Hackman et al., 2019). Tikhonov regularization prefers the solution with
the anti-solar differential rotation, since it tries to minimize the
temperature gradient. From Fig. 1 one can see that the model for the I
profiles is categorically deeper than the observations, which also indicates
that the I data is not modeled as well as the V data. Thus, we do not conclude
anything based on the differential rotation analysis of the I data, except
that the apparent results from I and V data can be radically different. A
$\sigma$ map for the I data is shown in Appendix A.

Figure 4: Magnetic field and brightness maps for all six nights with
differential rotation parameters $\alpha=0.075$ and $P_{\mathrm{rot}}=1.31$ d.
The observed phases are not indicated, since with differential rotation there
is no unique rotation period. The horizontal line indicates the limit of the
visible surface due to the inclination. The maps are shown for the epoch
HJD=2455699.0 Figure 5: The Stokes profiles for the best value for
differential rotation, $\alpha=0.075$ and $P_{\mathrm{rot}}=1.31$ d. The black
lines are the observed LSD profiles, and the red lines are fits made to these.
The rotational phase is indicated as HJD relative to 2455699.0 Figure 6: The
deviation $\sigma$ as a function of $\alpha$, where for each $\alpha$ the
value for $P_{\mathrm{rot}}$ which minimizes $\sigma$ is used.

Figure 7: Magnetic field and brightness maps with differential rotation parameters $\alpha=0.075$ and $P_{\mathrm{rot}}=1.31$ d, for subsets with nights 1-3 (up), nights 2-5 (centre) and nights 4-6 (down). The observed phases are not indicated, since with differential rotation the corresponding longitude is not trivial. The horizontal line indicates the limit of the visible surface due to the inclination. All maps are shown for the epoch HJD=2455699.0 Table 3: Deviation between the model and observations $\sigma$ for different data sets and different differential rotation parameters for Stokes V data. Nights included | $\alpha$ | $P_{\mathrm{rot}}$ | $\sigma\times 10^{-5}$
---|---|---|---
1-6 | 0 | 1.3300 | 3.3560
1-3 | 0 | 1.3300 | 1.9456
2-5 | 0 | 1.3300 | 2.1631
4-6 | 0 | 1.3300 | 2.8015
1-6 | 0.075 | 1.31 | 2.5298
1-3 | 0.075 | 1.31 | 1.9440
2-5 | 0.075 | 1.31 | 2.0933
4-6 | 0.075 | 1.31 | 2.6936

###  6.3 Case 3: Abruptly evolving magnetic field and differential rotation

The two scenarios, abrupt evolution of the magnetic field and differential
rotation, do not exclude each other. We thus still explore the possibility of
implementing them both.

With the best fit value for the differential rotation, we tried to repeat the
analysis keeping the parameters $\alpha=0.075$ and $P_{\mathrm{rot}}=1.31$ d,
but now inserting them separately to the subsets 1-3, 2-5 and 4-6 of the data.
As is seen from Table 3, the inclusion of differential rotation does not much
affect $\sigma$ for subset 1-3, while subsets 2-5 and 4-6 now have a clearly
better fit. The ZDI maps of the three subsets, including differential
rotation, are shown in Fig. 7. The main features of the maps appear fairly
similar as in the maps without differential rotation. Now, though, the
magnetic field strength increases from the subset 1-3 to subsets 2-5 and 4-6,
while in the maps without differential rotation it decreased from subsets 1-3
and 2-5 to subset 4-6. The corresponding Stokes LSD profiles are shown in
Appendix B.

##  7 Discussion

It is evident that a constant magnetic field without differential rotation
cannot explain our data. A comparison of the two different alternatives, a
rapidly changing magnetic field or solar-like differential rotation, presented
in Table 3, reveals that the values for the deviation between model and
observations, $\sigma$, do not clearly prefer one alternative over the other.
For the subset with nights 1-3, $\sigma$ is very similar with or without
differential rotation, whereas for subset 4-6 the solution with differential
rotation is preferred. It seems difficult to clearly distinguish between the
two cases. The fit to the LSD profiles without magnetic field evolution (Fig.
5) is fairly good, and the inclusion of differential rotation does also
decrease $\sigma$ for the models with magnetic field evolution (although only
marginally for nights 1-3). Thus it seems that differential rotation is
required to explain the observations. We would argue that the most likely
explanation is a situation were differential rotation is present, and some
abrupt evolution does indeed also happen in the magnetic field configuration.
Indeed, if nights 1-3 and 4-6 are modeled independently, including
differential rotation, the resulting magnetic field strengths are different,
which is the biggest evidence for the abrupt magnetic field evolution to be
present besides the differential rotation. At longitude 120, there is a radial
magnetic field feature in the 4-6 map and possibly 2-5 map, which is missing
in the 1-3 map. If this feature is real, its appearance could also explain the
observed difference.

The deviation $\sigma$ increases gradually in all models from nights 1-3 to
2-5 and 4-6. The appearing magnetic region could be an explanation for this,
but it could also be explained with the S/N, which is lower for the last
nights, especially night 6, and will thus introduce larger residuals between
the model and observations. Thus we cannot say with certainty, how significant
the possible magnetic field evolution would be.

##  8 Implications for ZDI

In this section we assume that the hypothesis of the abruptly evolving
magnetic field is true, and discuss its possible implications.

Since this data set is of unusually good quality, when compared to other sets
collected for V889 Her, it is impossible to say if it is common for the
magnetic field of V889 Her to evolve this rapidly. At least this result
demonstrates that (Z)DI maps should be regarded with some caution, and the
underlying assumption that the magnetic field stays constant through the time
span of the observations should be kept in mind. For example, the observations
in Willamo et al. (2019) were generally done over a period of 11 nights, an
almost two times longer time span than the six nights analyzed here.
Differential rotation was not taken into account in that study either. Its
main result, a persistent polar spot, should not be too much affected by
abrupt magnetic field evolution or differential rotation, though.

According to the photometry analyzed in Willamo et al. (2019), V889 Her was in
a phase of increasing activity in 2011. Around this time, the fairly regular
cycle was replaced with a monotonous increase of activity, lasting at least
until 2019. It would be reasonable to expect that this kind of abrupt changes
in the magnetic field are more common in times of higher activity, such as for
the Doppler images published by Willamo et al. (2019) during 2012-2017.

The possibility of rapid magnetic field evolution should be kept in mind for
other stars as well. Still, it would be reasonable to assume that the magnetic
field of V889 Her might change considerably rapidly, since V889 Her is a very
active ($\log R^{\prime}_{\mathrm{HK}}=-4.175$; Lehtinen et al. (2016)),
rapidly rotating star, with a rotation period of approximately 1.3 days. Also
in the study by Lehtinen et al. (2022) of LQ Hya, a similar, almost as rapidly
rotating star ($P_{\mathrm{rot}}\approx 1.6$ d), there are a few spectra at
close by rotational phases (separated by some stellar rotations), which
display some variation in the Stokes V profiles, although much less than in
the data analyzed here. In Donati et al. (2016, 2017), ZDI was applied to the
T Tauri star V830 Tau ($P_{\mathrm{rot}}\approx 2.7$ d) using two
approximately one month long data sets, separated by one month. In this case,
there were slight variations in the magnetic topology, which were assigned to
differential rotation, but essentially a stable magnetic surface structure and
no apparent problems with ZDI, even for such a long observation span. With
less active, more slowly rotating stars, these kind of rapid variations might
be less common. Nevertheless, this suggests that rapid variations in the
magnetic fields of active stars can happen, and more caution should be applied
in the interpretation of results from ZDI. Especially if differential rotation
is not taken into account, both of these possible effects can affect the
results.

##  9 Field dynamics

Table 4: The amounts of poloidal and axisymmetric magnetic field, and the average and maximum magnetic field strength for each model. Multiple definitions for the axisymmetry are shown (with modes $m<\ell/2$ and $m=0$ being counted as axisymmetric). The amounts of toroidal and non-axisymmetric magnetic field can be calculated from these as 100% - the poloidal or axisymmetric component. Nights included | $\alpha$ | $P_{\mathrm{rot}}$ |  $E_{\mathrm{pol}}$ [%] |  $E_{m<\ell/2}$ [%] |  $E_{m=0}$ [%] |  $\langle|B|\rangle$ [G] |  $B_{\mathrm{max}}$ [G]
---|---|---|---|---|---|---|---
1-6 | 0 | 1.3300 | 94 | 93 | 90 | 1050 | 2720
1-3 | 0 | 1.3300 | 92 | 89 | 86 | 1160 | 3310
2-5 | 0 | 1.3300 | 96 | 98 | 96 | 1450 | 3430
4-6 | 0 | 1.3300 | 80 | 75 | 59 | 540 | 2570
1-6 | 0.075 | 1.31 | 98 | 98 | 97 | 2030 | 4730
1-3 | 0.075 | 1.31 | 92 | 90 | 86 | 1030 | 3350
2-5 | 0.075 | 1.31 | 98 | 98 | 97 | 1910 | 4450
4-6 | 0.075 | 1.31 | 93 | 86 | 83 | 1780 | 4950

Regardless of the model for differential rotation or magnetic field evolution,
the magnetic field is very much dominated by its poloidal and axisymmetric
components. These are listed in Table 4, along with the average and maximum
field strength. The magnetic field is notably stronger (average field strength
between 1-2 kG in most scenarios444The magnetic field strength might exceed
the limit for the weak field approximation ($\approx$ 1 kG) that is assumed
for the LSD profiles (see Kochukhov et al., 2010), which means that the field
strength in regions with the strongest field is not reliable. The order of
magnitude should, nevertheless, be correct.) than in the data from 2004 and
2005 from Marsden et al. (2006) and Jeffers and Donati (2008), where the field
strength was much below 1 kG. Those magnetic field maps also have a much
stronger azimuthal component than our maps. Some changes in the magnetic field
topology may thus have occurred. During 2004-2005 V889 Her was close to a spot
activity maximum, according to the photometric cycle seen in Willamo et al.
(2019), but in 2011 the star had entered a state of higher activity than the
previous maximum, as its cyclicity had temporarily been replaced with a
monotonous decrease in brightness. The state of higher spot activity may thus
explain the stronger magnetic field and different magnetic topology.

It should be noted, that the very axisymmetric magnetic field topology makes
the estimation of differential rotation more difficult, and could lead to
uncertainties regarding it. Differential rotation is also otherwise known to
be a difficult parameter to estimate reliably.

The brightness map, again in all the models for differential rotation or
magnetic field evolution, is dominated by a high-latitude spot, which is
centered a bit off the pole. The low-latitude features are different in the
different models, as is the contrast of the high-latitude spot, but the main
result from the brightness map, the persistent polar or high-latitude spot, is
the same as in earlier Doppler imaging studies of V889 Her (e.g. Willamo et
al., 2019).

Assuming some correlation between spots and magnetic fields, the possible
rapid field evolution would indicate that there are spots appearing and
disappearing rapidly on V889 Her. Giles et al. (2017) found on their sample of
Kepler stars, that hotter stars have shorter spot life times than cooler
stars; in their Figure 8, there appear to be spots with short life times, with
decay times around 10 days (which they also used as a lower limit for the
decay time) only in G and F-type stars, while not in the cooler K and M-stars.
As a G-type star, the possible rapid field evolution of V889 Her would fit in
this picture.

The most dramatic topological changes, seen in the radial and meridional
magnetic field, happen at low latitudes. Note that these changes still take
place above the equator – with the spherical harmonic decomposition of the
magnetic field, where the sum of the magnetic field penetrating the stellar
surface is forced to be zero, the magnetic features below the equator are
quite uncertain. This is because we do not have any information of the
magnetic field around the invisible pole, and the inversion forces a strong
lower latitude magnetic field with opposite polarity to compensate for the
observed magnetic field around the visible pole, even though a large part of
the magnetic field is probably unobservable to us.

If the possible magnetic field evolution truly happens close to equatorial
regions, as the radial feature around longitude 120 could indicate, then one
phenomenon, seen in numerical magnetohydrodynamic simulations, which could
possibly give an explanation for this, are ‘magnetic wreaths’. They are
banded, large-scale magnetic structures, which arise at low latitudes on both
hemispheres in simulations of rapidly rotating stars (Brown et al., 2010).
More recently, higher-resolution simulations (see, e.g., Viviani et al., 2018)
have revealed that these structures can exhibit small-scale structures, as the
toroidal band of magnetic field is twisted and obtains a complicated helical
structure. This means that it is not only manifested in the azimuthal magnetic
field component, but also in the radial and meridional ones. The wreaths show
faster and more dynamical behavior than the non-axisymmetric large-scale
structures appearing at higher latitudes in these simulations. Since it is
known that there is a large contribution of small-scale fields, which is not
captured by ZDI (Kochukhov et al., 2020), the rapid time scale of the magnetic
field evolution supports the conclusion that the actual changes in the
magnetic field happen on small scales. Since these wreaths have strong
structures of small-scale fields, they could be expected to evolve on short
time scales, due to rapid processes such as reconnection and diffusion. Thus,
they provide a possible explanation for the fast low latitude changes seen in
V889 Her.

##  10 Conclusions

We have applied ZDI for an extensive set of spectropolarimetry from six
consecutive nights in May 2011 for V889 Her. We showed that a constant
magnetic field on a solidly rotating star cannot explain the data. Thus we
studied two alternative scenarios; rapid evolution of the magnetic field and
differential rotation. From our analysis we conclude that it is not easy to
distinguish between these two alternatives from ZDI alone. It is evident that
strong differential rotation is present in V889 Her. The abrupt evolution of
the magnetic field is not as certain, but it seems likely that at least some
evolution has also occurred. Our result disagrees with the theory and models
where rapidly rotating stars do not show strong differential rotation (compare
our $\alpha=0.075$ to e.g. $\alpha<0.01$ for simulations with comparable
rotation (Viviani et al., 2018), or $\alpha\lesssim 0.05$ for Kepler stars
with $T_{\mathrm{eff}}<6000$ K in Reinhold and Gizon (2015)). We would also
urge to caution, when interpreting results of ZDI from long data sets, since
it seems possible that some evolution of the large-scale magnetic fields of
rapid rotators is possible in short time-scales of only a few days.

###### Acknowledgements.

Based on observations made with the HARPSpol instrument on the ESO 3.6 m
telescopes at the La Silla Observatory under programme ID 087.D-0771(A). T.W.
acknowledges financial support of the Väisälä and von Frenckell foundations.
T.H. acknowledges the financial support from the Academy of Finland for the
project SOLSTICE (decision No. 324161). M.J.K-L. acknowledges the support of
the Academy of Finland ReSoLVE Centre of Excellence (grant number 307411).
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(Project UniSDyn, grant agreement n:o 818665). O.K. acknowledges support by
the Swedish Research Council and the Swedish National Space Agency. We thank
the anonymous referee for comments that helped to improve the paper.

##  Appendix A Differential rotation for Stokes I data

Figure 8: Map of the deviation $\sigma$ for differential rotation parameters
of V889 Her, using Stokes I data. The colour scaling corresponds to the
deviation between the model and observations in the respective ZDI map, with
the unit being $10^{-4}$. The best fit is indicated with a cross.

##  Appendix B Stokes profiles for the model with both differential rotation
and magnetic field evolution

Figure 9: The Stokes profiles for the best value for differential rotation,
$\alpha=0.075$ and $P_{\mathrm{rot}}=1.31$ d, when the magnetic field is also
allowed to evolve. These are independently modeled for nights 1-3 (left),
nights 2-5 (centre) and nights 4-6 (right). The black lines are the observed
Stokes LSD profiles, and the red lines are fits made to these.

## References

  * Brown et al. (2010) Benjamin P. Brown, Matthew K. Browning, Allan Sacha Brun, Mark S. Miesch, and Juri Toomre.  Persistent Magnetic Wreaths in a Rapidly Rotating Sun.  _ApJ_ , 711(1):424–438, March 2010.  10.1088/0004-637X/711/1/424.
  * Brown et al. (1991) S. F. Brown, J. F. Donati, D. E. Rees, and M. Semel.  Zeeman-Doppler imaging of solar-type and AP stars. IV. Maximum entropy reconstruction of 2D magnetic topologies.  _A &A_, 250:463, October 1991.
  * Donati et al. (2016) J. F. Donati, C. Moutou, L. Malo, C. Baruteau, L. Yu, E. Hébrard, G. Hussain, S. Alencar, F. Ménard, J. Bouvier, P. Petit, M. Takami, R. Doyon, and A. Collier Cameron.  A hot Jupiter orbiting a 2-million-year-old solar-mass T Tauri star.  _Nature_ , 534(7609):662–666, June 2016.  10.1038/nature18305.
  * Donati et al. (2017) J. F. Donati, L. Yu, C. Moutou, A. C. Cameron, L. Malo, K. Grankin, E. Hébrard, G. A. J. Hussain, A. A. Vidotto, S. H. P. Alencar, R. D. Haywood, J. Bouvier, P. Petit, M. Takami, G. J. Herczeg, S. G. Gregory, M. M. Jardine, J. Morin, and MaTYSSE Collaboration.  The hot Jupiter of the magnetically active weak-line T Tauri star V830 Tau.  _MNRAS_ , 465(3):3343–3360, March 2017.  10.1093/mnras/stw2904.
  * Frasca et al. (2010) A. Frasca, K. Biazzo, Z. Kővári, E. Marilli, and Ö. Çakırlı.  Photospheric and chromospheric activity on the young solar-type star HD 171488 (V889 Herculis).  _A &A_, 518:A48, July 2010.  10.1051/0004-6361/201014460.
  * Giles et al. (2017) Helen A. C. Giles, Andrew Collier Cameron, and Raphaëlle D. Haywood.  A Kepler study of starspot lifetimes with respect to light-curve amplitude and spectral type.  _MNRAS_ , 472(2):1618–1627, December 2017.  10.1093/mnras/stx1931.
  * Hackman et al. (2016) T. Hackman, J. Lehtinen, L. Rosén, O. Kochukhov, and M. J. Käpylä.  Zeeman-Doppler imaging of active young solar-type stars.  _A &A_, 587:A28, March 2016.  10.1051/0004-6361/201527320.
  * Hackman et al. (2019) T. Hackman, I. Ilyin, J. J. Lehtinen, O. Kochukhov, M. J. Käpylä, N. Piskunov, and T. Willamo.  Starspot activity of HD 199178. Doppler images from 1994-2017.  _A &A_, 625:A79, May 2019.  10.1051/0004-6361/201834763.
  * Huber et al. (2009) K. F. Huber, U. Wolter, S. Czesla, J. H. M. M. Schmitt, M. Esposito, I. Ilyin, and J. N. González-Pérez.  Long-term stability of spotted regions and the activity-induced Rossiter-McLaughlin effect on V889 Herculis. A synergy of photometry, radial velocity measurements, and Doppler imaging.  _A &A_, 501:715–728, July 2009.  10.1051/0004-6361/200810638.
  * Järvinen et al. (2008) S. P. Järvinen, H. Korhonen, S. V. Berdyugina, I. Ilyin, K. G. Strassmeier, M. Weber, I. Savanov, and I. Tuominen.  Magnetic activity on V889 Herculis. Combining photometry and spectroscopy.  _A &A_, 488:1047–1055, September 2008.  10.1051/0004-6361:200809837.
  * Jeffers and Donati (2008) S. V. Jeffers and J.-F. Donati.  High levels of surface differential rotation on the young G0 dwarf HD171488.  _MNRAS_ , 390:635–644, October 2008.  10.1111/j.1365-2966.2008.13695.x.
  * Jetsu and Pelt (1999) L. Jetsu and J. Pelt.  Three stage period analysis and complementary methods.  _A &AS_, 139:629–643, November 1999.  10.1051/aas:1999411.
  * Kővári et al. (2011) Z. Kővári, A. Frasca, K. Biazzo, K. Vida, E. Marilli, and Ö. Çakırlı.  Differential rotation on the young solar analogue V889 Herculis.  In D. Prasad Choudhary and K. G. Strassmeier, editors, _Physics of Sun and Star Spots_ , volume 273 of _IAU Symposium_ , pages 121–125, August 2011.  10.1017/S1743921311015110.
  * Kitchatinov and Rüdiger (1999) L. L. Kitchatinov and G. Rüdiger.  Differential rotation models for late-type dwarfs and giants.  _A &A_, 344:911–917, April 1999.
  * Kochukhov et al. (2010) O. Kochukhov, V. Makaganiuk, and N. Piskunov.  Least-squares deconvolution of the stellar intensity and polarization spectra.  _A &A_, 524:A5, December 2010.  10.1051/0004-6361/201015429.
  * Kochukhov et al. (2014) O. Kochukhov, T. Lüftinger, C. Neiner, E. Alecian, and MiMeS Collaboration.  Magnetic field topology of the unique chemically peculiar star CU Virginis.  _A &A_, 565:A83, May 2014.  10.1051/0004-6361/201423472.
  * Kochukhov et al. (2020) O. Kochukhov, T. Hackman, J. J. Lehtinen, and A. Wehrhahn.  Hidden magnetic fields of young suns.  _A &A_, 635:A142, March 2020.  10.1051/0004-6361/201937185.
  * Kochukhov (2016) Oleg Kochukhov.  _Doppler and Zeeman Doppler Imaging of Stars_ , volume 914, page 177.  2016\.  10.1007/978-3-319-24151-7_9.
  * Kupka et al. (1999) F. Kupka, N. Piskunov, T. A. Ryabchikova, H. C. Stempels, and W. W. Weiss.  VALD-2: Progress of the Vienna Atomic Line Data Base.  _A &AS_, 138:119–133, July 1999.  10.1051/aas:1999267.
  * Lehtinen et al. (2016) J. Lehtinen, L. Jetsu, T. Hackman, P. Kajatkari, and G. W. Henry.  Activity trends in young solar-type stars.  _A &A_, 588:A38, April 2016.  10.1051/0004-6361/201527420.
  * Lehtinen et al. (2022) J. J. Lehtinen, M. J. Käpylä, T. Hackman, O. Kochukhov, T. Willamo, S. C. Marsden, S. V. Jeffers, G. W. Henry, and L. Jetsu.  Topological changes in the magnetic field of LQ Hya during an activity minimum.  _A &A_, 660:A141, April 2022.  10.1051/0004-6361/201936780.
  * Marsden et al. (2006) S. C. Marsden, J.-F. Donati, M. Semel, P. Petit, and B. D. Carter.  Surface differential rotation and photospheric magnetic field of the young solar-type star HD 171488 (V889 Her).  _MNRAS_ , 370:468–476, July 2006.  10.1111/j.1365-2966.2006.10503.x.
  * Montes et al. (2001) D. Montes, J. López-Santiago, M. C. Gálvez, M. J. Fernández-Figueroa, E. De Castro, and M. Cornide.  Late-type members of young stellar kinematic groups - I. Single stars.  _MNRAS_ , 328:45–63, November 2001.  10.1046/j.1365-8711.2001.04781.x.
  * Piskunov et al. (2011) N. Piskunov, F. Snik, A. Dolgopolov, O. Kochukhov, M. Rodenhuis, J. Valenti, S. Jeffers, V. Makaganiuk, C. Johns-Krull, E. Stempels, and C. Keller.  HARPSpol — The New Polarimetric Mode for HARPS.  _The Messenger_ , 143:7–10, March 2011.
  * Piskunov and Valenti (2002) N. E. Piskunov and J. A. Valenti.  New algorithms for reducing cross-dispersed echelle spectra.  _A &A_, 385:1095–1106, April 2002.  10.1051/0004-6361:20020175.
  * Piskunov et al. (1995) N. E. Piskunov, F. Kupka, T. A. Ryabchikova, W. W. Weiss, and C. S. Jeffery.  VALD: The Vienna Atomic Line Data Base.  _A &AS_, 112:525, September 1995.
  * Piskunov et al. (1990) N.E. Piskunov, I. Tuominen, and O. Vilhu.  Surface imaging of late-type stars.  _A &A_, 230:363–370, 1990.
  * Reiners and Schmitt (2003) A. Reiners and J. H. M. M. Schmitt.  Differential rotation in rapidly rotating F-stars.  _A &A_, 412:813–819, December 2003.  10.1051/0004-6361:20034255.
  * Reinhold and Gizon (2015) Timo Reinhold and Laurent Gizon.  Rotation, differential rotation, and gyrochronology of active Kepler stars.  _A &A_, 583:A65, November 2015.  10.1051/0004-6361/201526216.
  * Rosén et al. (2016) L. Rosén, O. Kochukhov, T. Hackman, and J. Lehtinen.  Magnetic fields of young solar twins.  _A &A_, 593:A35, September 2016.  10.1051/0004-6361/201628443.
  * Strassmeier and Bartus (2000) K. G. Strassmeier and J. Bartus.  Doppler imaging of stellar surface structure. XII. Rapid spot changes on the RS CVn binary V711 Tauri = HR 1099.  _A &A_, 354:537–550, February 2000.
  * Strassmeier et al. (2003a) K. G. Strassmeier, L. Kratzwald, and M. Weber.  Doppler imaging of stellar surface structure. XX. The rapidly-rotating single K2-giant HD 31993 = V1192 Orionis.  _A &A_, 408:1103–1113, September 2003a.  10.1051/0004-6361:20031029.
  * Strassmeier et al. (2003b) K. G. Strassmeier, T. Pichler, M. Weber, and T. Granzer.  Doppler imaging of stellar surface structure. XXI. The rapidly-rotating solar-type star HD 171488 = V889 Hercules.  _A &A_, 411:595–604, December 2003b.  10.1051/0004-6361:20031538.
  * Viviani et al. (2018) M. Viviani, J. Warnecke, M. J. Käpylä, P. J. Käpylä, N. Olspert, E. M. Cole-Kodikara, J. J. Lehtinen, and A. Brandenburg.  Transition from axi- to nonaxisymmetric dynamo modes in spherical convection models of solar-like stars.  _A &A_, 616:A160, August 2018.  10.1051/0004-6361/201732191.
  * Vogt et al. (1987) Steven S. Vogt, G. Donald Penrod, and Artie P. Hatzes.  Doppler images of rotating stars using maximum entropy image reconstruction.  _ApJ_ , 321:496–515, 1987.
  * Willamo et al. (2019) T. Willamo, T. Hackman, J. J. Lehtinen, M. J. Käpylä, I. Ilyin, G. W. Henry, L. Jetsu, O. Kochukhov, and N. Piskunov.  Long-term spot monitoring of the young solar analogue V889 Herculis.  _A &A_, 622:A170, February 2019.  10.1051/0004-6361/201834562.
  * Willamo et al. (2022) T. Willamo, J. J. Lehtinen, T. Hackman, M. J. Käpylä, O. Kochukhov, S. V. Jeffers, H. Korhonen, and S. C. Marsden.  Zeeman-Doppler imaging of five young solar-type stars.  _A &A_, 659:A71, March 2022.  10.1051/0004-6361/202141649.