f-mode strengthening from a localised bipolar subsurface magnetic field

ABSTRACT Recent numerical work in helioseismology has shown that a periodically varying subsurface magnetic field leads to a fanning of the f-mode, which emerges from a density jump at the surface. In an attempt to model a more realistic situation, we now modulate this periodic variation with an envelope, giving thus more emphasis on localised bipolar magnetic structures in the middle of the domain. Some notable findings are: (i) compared to the purely hydrodynamic case, the strength of the f-mode is significantly larger at high horizontal wavenumbers k, but the fanning is weaker for the localised subsurface magnetic field concentrations investigated here than the periodic ones studied earlier; (ii) when the strength of the magnetic field is enhanced at a fixed depth below the surface, the fanning of the f-mode in the diagram increases proportionally in such a way that the normalised f-mode strengths remain nearly the same in different such cases; (iii) the unstable Bloch modes reported previously in case of harmonically varying magnetic fields are now completely absent when more realistic localised magnetic field concentrations are imposed beneath the surface, thus suggesting that the Bloch modes are unlikely to be supported during most phases of the solar cycle; (iv) the f-mode strength appears to depend also on the depth of magnetic field concentrations such that it shows a relative decrement when the maximum of the magnetic field is moved to a deeper layer. We argue that detections of f-mode perturbations such as those being explored here could be effective tracers of solar magnetic fields below the photosphere before these are directly detectable as visible manifestations in terms of active regions or sunspots.


Introduction
The Sun supports a wide variety of waves that carry useful information about the internal solar structure, which can be inferred by employing the methods of helioseismology (see, e.g. Gough 1987;Christensen-Dalsgaard 2003). Local analysis using especially the surface gravity mode (or f -mode) is useful in studying the near-surface structure (e.g. Hanasoge et al. 2008;Daiffallah et al. 2011;Felipe et al. 2012Felipe et al. , 2013. Some properties of surface waves in idealised settings involving magnetic fields were explored in Chandrasekhar (1961), Roberts (1981), Roberts (1989, 1992) and Miles et al. (1992). It is of great interest to probe interior magnetic fields of the Sun using the techniques of helioseismology (see Thompson 2006, for a review on the subject of magnetohelioseismology). Systematic changes, notably in the frequencies of the global f -mode, as a function of solar cycle were found and discussed in detail in Thompson (2006) and Pintér (2008). Thompson (2006) found evidence of a 500 G magnetic field at a depth of about 5 Mm and suggested 2% modulation in turbulent convection velocities from solar minimum to maximum, based on the observed cycle-dependence of mean frequency shifts of the f and p-modes. Moreover, it is reasonable to expect magnetically induced variations in the f -mode on much shorter timescales, e.g. during localised magnetic flux emergence leading to the formation of active regions (ARs) or sunspots.
Much of the earlier studies on the global f -mode of the Sun focussed primarily on the frequency shifts that were observed (Libbrecht et al. 1990;Fernandes et al. 1992). The frequencies were significantly smaller than the theoretically expected values, where both the shift and line width grow with the spherical harmonic degree. Subsequent studies explained these findings by invoking turbulent background motions and deriving a generalised dispersion relation of the f -mode in the presence of a random velocity field Roberts 1993a, 1993b;Medrek et al. 1999;Murawski 2000aMurawski , 2000bMole et al. 2008). The influence of coherent (for modelling supergranulation) as well as random (mimicking near-surface granulation) flows on the f -mode were explored in detail by Murawski (2000a) where it was found that, while a space-dependent random flow causes a decrement, a time-dependent random flow can enhance the frequencies. Observations were thus explained in terms of the parameters of the chosen velocity field. However, these studies ignored the effect of magnetic fields, which can increase the f -mode frequencies (e.g. Chandrasekhar 1961;Roberts 1981;Miles and Roberts 1992;Miles et al. 1992).
In a series of works, Bogdan (1993, 1997) and Cally et al. (1994) investigated the interaction of f and p-modes with a vertical magnetic field and found that a partial conversion of these modes into slow magnetoacoustic modes takes place whenever they encounter a vertical field resembling those of sunspots. Parchevsky and Kosovichev (2009) numerically explored the effects of inclined magnetic fields and noted that the f -modes are more strongly affected by the background magnetic field than the p-modes. Singh et al. (2015) studied numerically various properties of the f as well as p and g-modes in a wide variety of magnetic backgrounds. They found that horizontal magnetic fields cause an increase in the f -mode frequencies -as expected. But their dependencies are more complicated in the presence of vertical or oblique magnetic fields, which may be more relevant for the predominantly vertical fields of sunspots. In this case, the f -mode frequencies are enhanced relative to their nonmagnetic values at intermediate horizontal wavenumbers, but decreased at large wavenumbers.
In the presence of magnetic fields in the solar atmosphere, Alfvén and magnetosonic waves are known to couple resonantly with the global oscillations, affecting mainly the frequencies, line widths, and penetration of the f and p-modes into the solar atmosphere (see Erdélyi 2006, for a review and references therein). Another important finding was that resonant interaction between global modes and Alfvén waves causes a damping of the f and p-modes due to dissipative effects near the resonance frequency. Pintér and Erdélyi (2018) further study this situation by deriving the dispersion relation of the f -mode in a magnetically coupled solar interior-atmosphere system to obtain its frequency shifts. For a magnetised atmosphere, these were found to be positive relative to an unmagnetised one.
It is of great interest to use numerical simulations to study the effects of subsurface magnetic fields on both the acoustic or p-modes and the f -mode. Such studies aim at refining the interpretation of helioseismic measurements which use sound waves to infer the internal structure of the Sun and its internal motions (e.g. Basu 2016;Hanasoge et al. 2016). This is necessary, because magnetic fields complicate the usage of helioseismic inversion techniques, as their presence gives rise to the modification of sound and gravity waves into magnetoacoustic and magnetogravity ones which are difficult to account for (Thomas, 1983;Campos, 2011). Furthermore, if the properties and behaviour of p and f -modes in the presence of magnetic fields of different strength, topology or location, in particular depth, were known well enough from the numerical models, we might be able to infer the subsurface magnetic fields from helioseismic measurements. From such inferences, we may be able to learn about the origin of subsurface magnetic fields, that is, about the solar dynamo mechanism (Brandenburg 2005;Charbonneau 2010;Käpylä et al. 2012), and about the process of concentrating magnetic fields into sunspots and ARs (see, e.g. Brandenburg et al. 2016;Käpylä et al. 2016 for a mechanism enabling a self-consistent formation of localised magnetic flux concentrations).
Such numerical simulations, enabling us to study the effects of magnetic fields on the naturally occurring modes of oscillations, have recently been performed with the Pencil Code 1 . Singh et al. (2015) introduced a modelling framework with a piecewise isothermal atmosphere, where the upper layer is mimicking a hot corona and the lower one a (cooler) convection zone. The two layers are separated by a jump in density and temperature, which represents the solar surface and enables the presence of the f -mode. For a simple approximation of the convective turbulence, random hydrodynamic forcing was applied in the lower layer. In such a setup, acoustic (p), internal gravity (g), and surface gravity (f ) modes are all self-consistently driven, in contrast to other types of approaches, which selectively produce the modes based on linearised equations (Daiffallah et al. 2011;Schunker et al. 2011). Singh et al. (2015 studied the influence of uniformly imposed magnetic fields on the p and f -mode properties, and verified the expectation of f -modes being sensitive to the presence of magnetic fields and p-modes being less affected. The major effect on the f -modes was an increase in their mode frequencies. Singh et al. (2014) developed the setup further to include non-uniform magnetic fields with harmonic profiles within the convection zone. This work revealed the fanning effect of the f -mode, that is, an increase of its line width towards higher wavenumbers. The width of the fan and its asymmetry could be directly related to the strength and location of the magnetic field.
These numerical studies and their predictions led to an observational case study of f -modes in relation to the emergence of about half a dozen ARs, observed with HMI, and the resulting line-of-sight Dopplergrams and magnetograms (Singh et al. 2016). This study reported strengthening of the local f -mode about two days before the emergence of ARs at the same corotating patch. It was noted that such a precursor signal can be detected by isolating high-degree f -modes through careful fitting and subsequent subtraction of background and p-modes. It was argued there that the precursor signal is best seen when an AR forms in isolation, i.e. far from other existing ARs which can "pollute" the signal. It is known that the sunspots or ARs absorb the f -mode power, thus causing its damping (Cally and Bogdan 1997), and therefore it is expected to be harder to extract the precursor signal associated with a newly forming AR in a "crowded" environment with many existing ARs. A possible cleaning procedure to still extract the signal was discussed in Singh et al. (2016). Although this is yet to be confirmed for larger data sets and with independent observational techniques, the potential significance of such f -mode related precursors cannot be understated. Based on the photospheric velocity measurements, Khlystova and Toriumi (2017) detected plasma upflows somewhat before the emergence of two small ARs. A number of previous case studies have reported detections of subphotospheric velocities associated with emerging ARs using different techniques, such as, helioseismic holography, ring-diagram or time-distance analysis (Komm et al. 2008;Hartlep et al. 2011;Ilonidis et al. 2011;Birch et al. 2013; Barnes et al. 2014). Given that the f -mode eigenfunctions extend to depths of about a few Mm, they might be expected to be sensitive to such velocity perturbations.
In this paper we extend the model of Singh et al. (2014Singh et al. ( , 2015 to study the f -mode strengthening and fanning using inhomogeneous magnetic fields of different strength, topology, and location in the convection zone. In addition to using harmonic profiles which stretch over the whole horizontal extent, we use localised harmonic perturbations, to better mimic isolated active regions. We describe our model and basic definitions in section 2 and present our results in section 3. Conclusions are given in section 4.

Model setup and analysis technique
We consider here a model that is similar to that studied in Singh et al. (2014Singh et al. ( , 2015. The p, g and f -modes are excited in a self-consistent manner in a two-dimensional Cartesian x-z domain with a piecewise isothermal medium consisting of a cool lower layer, called bulk, and a hotter upper layer, called corona. The thicknesses of these layers are L d and L u , respectively, where subscripts "d" and "u" refer to the layers downward and upward of the interface. A non-uniform magnetic background field is maintained by including an electromotive force (EMF) in the uncurled version of the induction equation for the magnetic vector potential A. We solve the basic hydromagnetic equations, Du

∂A ∂t
where ρ is the density, u is the velocity, D/Dt = ∂/∂t + u·∇ is the advective time derivative, f is a forcing function (as specified in Brandenburg 2001) to drive seismic modes of low Mach number 2 , g = (0, 0, −g) is the gravitational acceleration, S ij = 1 2 (u i,j + u j,i ) − 1 3 δ ij ∇·u is the traceless rate of strain tensor, with commas denoting partial differentiation, ν = const is the kinematic viscosity, B = ∇ × A is the magnetic field, J = μ −1 0 ∇ × B is the current density, μ 0 is the vacuum permeability, s is the specific entropy, γ = c p /c v is the ratio of specific heats at constant pressure and volume, T is the temperature, E 0 is an external EMF specified below, and η = const is the magnetic diffusivity. Since the forcing amplitude is small, no additional hyperdiffusion is needed, and also no diffusion in the energy equation was included. The fluid is assumed to obey the equation of state of an ideal gas, hence the pressure is given by p = (c p − c v )ρT = ρc 2 s /γ . The medium is vertically stratified under constant gravity, g > 0, where we identify x and z as horizontal and vertical directions, respectively. All calculations are performed with the Pencil Code.
The variations of background density, pressure scale height, and pressure as functions of z are shown in figure 1. Similar to earlier works (Singh et al. 2014(Singh et al. , 2015 we introduce a sharp jump in density at the interface z = 0 with ρ u (0) ρ d (0), along with corresponding jumps in temperature. The adiabatic sound speed c s is maintained by the last term in (3), which guarantees the relaxation to constant average temperatures T d and T u in either subdomain within a relaxation time τ c (constant throughout the domain). In this way, the interface is created and maintained and the f -mode is naturally enabled. It is also known as the free surface mode. As the density decreases exponentially with height z in an isothermally stratified medium, it is given by (The pressure and density scale heights are equal for an isothermal layer.) The sharp jump in the thermodynamic quantities is quantified by the ratio We enforce a steady magnetic field B 0 by applying a constant external EMF E 0 = E A,B 0y e y . Two different types of magnetic structures are considered in the magnetohydrodynamic (MHD) simulations, called models A and B. Model A is characterised by a harmonic variation in both spatial dimensions and model B by a localised sine wave where F TH (x; x 1 , x 2 ) is a smoothed top hat as a function of x, centred w.r.t. x. It is unity for x ∈ [x 1 , x 2 ] and smoothly goes to zero outside this interval; see, e.g. figures 4 and 5 for the profiles of the imposed magnetic background in its saturated state. The transitions at x 1 and x 2 are modelled with a third-order polynomial of width w ≈ 0.3H d . Note that the sustained magnetic field B 0 drives a large-scale flow which, in turn, acts back on the field, so that its shape is not simply determined by η∇ 2 B 0 + ∇ × E 0 = 0. The A models are very similar to those considered in Singh et al. (2014), whereas the B ones were tailored to model an emerging active region more realistically: given that ARs hardly ever show a periodic pattern in longitude, we restricted the horizontal extent of the subsurface magnetic structures in the B models, mimicking more closely the fields of bipolar regions. Both models allow us to explore the effects of subsurface magnetic fields on the f -mode, naturally occurring at the interface in our minimalistic setup. It is nevertheless sufficiently realistic for understanding the solar f -mode that is expected to be unaffected by the choice of the thermodynamic background state of the medium. The top and bottom boundaries were chosen to be stress free and perfectly conducting, whereas periodicity was assumed in the horizontal direction. The length scales and frequencies are normalised by L 0 = γ H d = c 2 sd /g and ω 0 = g/c sd , and the dimensionless variables are indicated by tildae, i.e.k x = k x L 0 ,ω = ω/ω 0 , and so on. From the vertical velocity u z at z = 0, we construct the diagnostic k x -ω diagram (kω diagram for short) by taking the Fourier transform of u z (x, 0, t), givingû z (k x , ω), which here has the same dimension as velocity. The kω diagrams, such as those shown in figure 2 are constructed from the dimensionless quantitỹ where the mass-weighted root-mean-squared (rms) velocity u 0 = ρ 2 u 2 d / ρ 2 is employed for normalisation and the angle brackets denote volume averaging over the bulk. We define the fluid Reynolds and Mach numbers as Re = u 0 /(νk f ) and Ma = u 0 /c sd , respectively, where we have chosenk f = 20 for the wavenumber of the hydrodynamic low amplitude nonhelical forcing in (2).
In order to characterise the strength of the f -mode, we determine the normalised mode strength μ f as a function ofk x by fitting a Lorentzian to the line profile of the f -mode along the frequency axis and subtracting the continuum. Here we define the mode strength as where P denotes the normalised excess amplitude of the f -mode over the continuum at k x . This yields a wavenumber-dependent mode strength in a manner similar to that used in Singh et al. (2016) where, however, a more involved fitting algorithm was used to isolate the f -mode, which lies much closer to the p-modes in that more realistic setting with an isentropic stratification. It is useful to define an integrated normalised mode amplitude of the f -mode as which is called here, for short, the relative mode amplitude.

Results
We present results from a suite of hydrodynamic (H) as well as MHD runs (Sets A and B) with different geometries of the imposed background magnetic field. We refer the reader to table 1 for the parameters of all simulations.

Hydrodynamic runs
We performed four different hydrodynamic runs, H1-H4, to assess the roles of viscosity and vertical extent of f . We first show the diagnostic kω diagram, which clearly reveals the p, g, and f -modes; see the top left panel of figure 2 for model H1. While the g-modes are confined to frequenciesω < 1, the p-modes lie above the lineω =k x , as already seen in Singh et al. (2015) and as expected in the piecewise isothermal setup being considered here. As the present paper is focussed on the f -mode and its interaction with the subsurface magnetic fields, we refer the reader to Singh et al. (2015) for more details on the properties of the p and g-modes. The f -mode in the hydrodynamic cases (referred to as the classical f -mode) appears close to the theoretical curve given by (see, e.g. Chandrasekhar 1961;Gough 1987). We determine the strength of the f -mode using the kω diagram as discussed in the previous section, and show the wavenumber dependence of its normalised mode strength μ f in the bottom left panel of figure 2. No error bars have been added in this and similar figures below, because the actual fitting errors are small. However, in view of the non-smooth nature of the curves, surely other imperfections are present like, e.g. the limited integration time of the runs. The mode strengths from models H1 and H2 are nearly the same, although the viscosity of H2 is twice as large as that of H1. This is also reflected in the relative mode amplitude f listed in table 1. In both these cases, the same hydrodynamic forcing was employed in the whole domain. Unlike in models H1 and H2, the f -mode strength decreases systematically with wavenumberk x when the forcing is restricted to layers below the interface, namely toz < 0 in model H3 but toz < −0.2 in model H4. Compared to H3, μ f is smaller in H4 at nearly all wavenumbers. In figure 2 we also show a snapshot of the vertical motions in a representative model (H3). This figure demonstrates that no large scale flow patterns emerge in the HD simulations.
We note that the mode damping can possibly occur due to a resonant coupling of atmospheric Alfvén and magnetosonic waves with the global oscillations of the Sun, in which case dissipative effects play a major role (see Pintér et al. 2007). This is not applicable to the present work as we study the effects of magnetic fields that are confined to regions below the surface with negligible leakage into the atmosphere above. Such cases involving magnetic fields are presented in the following subsections below. Table 1. Summary of all simulations.x-z domain: 8π × π; grid: 1024 × 320; Pr m = 1; q = 0.11; |f |/g = 10 −4 for all runs.k B x = 0.5,k B z = 2 and w z = 0.28 were chosen in (6) and (7); x B = (x 2 − x 1 )/L 0 .

MHD A-type runs
Here we investigate the effects of harmonically varying subsurface magnetic fields, maintained by the EMF E A 0y as given in (6), on the f -mode. We refer first to the kω diagrams for runs A1 and A4, shown in the top panels of figure 3. The vertical stripes at low frequencies correspond to the unstable Bloch modes, which arise due to the periodicity of the background magnetic field (see Singh et al. 2014). Comparing these kω diagrams with the one from the suitable hydrodynamic run H1 (see table 1) shown in figure 2, we see that the f -modes in the MHD runs are very different, thus already providing evidence that the f -mode is sensitive to subsurface magnetism in our setup. This, in a sense, confirms earlier analytical work (e.g. Roberts 1981), where larger frequencies of the f -mode are expected, although the fanning and associated mode strengthening were not anticipated then. Let us define the z-dependent rms Alfvén speed, v A (z), as where • x denotes averaging along the x direction. The vertical profiles of v A (z)/c s (z) are shown in the bottom right panel of figure 3 for all A runs. In the bottom left panel of figure 3, we show the wavenumber dependence of the normalised mode strength μ f of the f -mode for the A runs and the corresponding hydrodynamic run H1. All the input parameters for these runs are the same (see table 1), except for the magnetic field strength which is gradually increased from A1 to A4 while keeping its maximum at the same depth; see figure 3. Compared to the run H1, the mode strength is clearly larger in the A runs for allk x > 1.5. Interestingly, the mode strength is nearly the same for runs A1-A4, while the line width (or fanning) of the f -mode increases from A1 to A4, as can be inferred from the kω diagrams. Snapshots of the vertical components B z and u z are shown in figure 4 for models A1 and A4 in the saturated stage of the background magnetic fields, where their maxima lie more than one pressure scale height below the interface at z = 0. While no large scale flow is noticeable in the weakly magnetised case A1, we do see a systematic flow pattern in model A4 with a stronger magnetic concentration. Despite the appearance of systematic flows in these MHD runs, the normalised f -mode strength is found to be nearly the same in all the A runs, as can also be seen from the corresponding values of f (identical when rounded to three decimal places) listed in table 1.

MHD B-type runs
Next we turn to MHD B-type models, where localised magnetic concentrations are maintained below the surface by the EMF E B 0y , as given in (7). Their horizontal extent x B = x 2 − x 1 is thus restricted, with x B ≈ 4π (i.e. half the horizontal extent of the domain) for the two BI-type models BI1 and BI2 and x B ≈ 7π for model BII; see the bottom left panels of figures 5 and 6, and table 1 for details.
The kω diagrams for models BI2 and BII are shown in figure 7. Remarkably, the unstable Bloch modes seen as vertical stripes at low frequencies in the A-type runs are completely absent here. These were first found in Singh et al. (2014) and thought to arise due to the periodicity of the background magnetic field (Berton and Heyvaerts 1987). Given that the appearance of solar ARs during most phases of the magnetic activity cycle is hardly ever periodic in longitude, the Bloch modes are unlikely to be seen. Yet, during very active phases of the Sun, these modes are potentially observable in the kω or ring diagrams that are constructed from sufficiently large patches covering many ARs which show a regular pattern of alternating positive and negative vertical magnetic field along the azimuthal direction. On the other hand, given the drastic differences between our A-and B-type runs, we must acknowledge that the Bloch modes are rather sensitive to the exact matching of the wavelength and the size of the patch. This limits the observability of Bloch modes severely.
As for models A1-A4, a direct comparison of diagnostic diagrams shown in figures 2 and 7 reveals that the f -mode is significantly perturbed in the presence of subsurface  magnetic fields; compare, e.g. BI2 with the corresponding hydrodynamic case H3 where the hydrodynamic forcing is restricted toz < 0 (see also table 1). The f -mode in the MHD run BI2 appears to exist at all wavenumbers shown, whereas it is markedly damped at large wavenumbers in the corresponding hydrodynamic run H3. Moreover, the fanning of the f -mode, as discussed above (see also Singh et al. 2014), is much weaker here for localised magnetic field concentrations.
The wavenumber dependence of μ f for the two BI-runs and the relevant hydrodynamic run H3, as well as the vertical profiles of the normalised Alfvén speed for these MHD runs are shown in the top panels of figure 5. For a meaningful comparison, the f -mode strengths were determined from data sets with identical time-spans, which can, in general, affect the mode strengths. The label H3s, used in figure 5(a) signifies that the data are from run H3, but μ f is determined from a shorter time-series, matching the ones in BI1 and BI2.
From figure 5(a) we see that μ f for the MHD run BI1 with weak magnetic concentration is nearly the same as in case of run H3s. Note also that the maximum of the Alfvén-to-sound speed ratio in this case lies at a larger depth compared to the A runs; see figure 5(b) and table 1. This may be another reason why we do not find mode strengthening in BI1. However, when the magnetic field is increased and its maximum moves upwards to somewhat shallower depths as in case of BI2, we immediately find the f -mode strengthening at wavenumbersk x ≥ 1.8. Snapshots of B z and u z are shown in the bottom panels of figure 5 for model BI2. The horizontal oscillations in u z in the proximity of regions of strong magnetic field may seem to be numerical artifacts, but they do not look like the regular "ringing" one encounters when a simulation is underresolved. As they are resolved by a few mesh points per period, they may well be real. Comparing these with the relevant hydrodynamic run H3 (see figure 2) we find that the vertical motions are comparable and have no large scale feature at the interface. As expected, a systematic flow pattern is indeed seen in deeper layers around the maximum of the magnetic concentration, but, as was argued before based on the A runs, it is not clear whether such flows can affect the mode strength of the f -mode which is determined from vertical motions only at the interface.
In those of our experiments where the forcing is further restricted to layersz < −0.2, as in the MHD model BII or in the hydrodynamic run H4 (see table 1), we still find evidence of f -mode strengthening in the presence of subsurface magnetic fields, but at larger wavenumbers (in the case of BII fork x ≥ 2.5). This is demonstrated in figure 6, where snapshots of B z and u z are also shown. In both models, μ f first decreases withk x , and, while it saturates at a constant level beyond about 2.5 in case of BII, it shows a monotonic decrement in H4; see figure 6(a). Interestingly, the f -mode strengthening is thus discernible in model BII, but is not present in model BI1, although both have similar magnetic field concentrations at nearly the same depth; see table 1, and compare panel (b) of figures 5 and 6. This is likely due to the larger horizontal extent of the magnetic structure in the case of BII as compared to BI1.
Understanding the exact physical mechanism, responsible for the strengthening of the f -mode observed here, is beyond the scope of the present paper. Both the imposed magnetic field and the systematic flow driven by it, can, in principle, influence frequency and strength of the f -mode. For assessing the relative importance of these two factors in our studies, let us compare Runs H1 and A1. As can be seen from the upper right panel of figure 4, the systematic flow can hardly be seen near the interface, mostly because of the magnetic field being weak in case A1. So we might conclude that the considerable increase in the relative mode amplitude f is mainly due to the magnetic field. By contrast, Run A4 (figure 4, lower right panel) exhibits a systematic flow pattern dominating over the random component of the flow. Hence, it cannot be decided to what extent the systematic flow is affecting f . In Runs A1-A4, f is nearly the same, although the magnetic field is changing by a factor of two. As the magnetic field alone seems to strengthen the f -mode, this could be an indication of a counteracting effect of the systematic flow, which grows with the magnetic field, too. More work is needed to understand the role of a systematic or mean flow component on the strength of the f -mode. This will be the focus of a separate investigation.

Conclusions
We have shown numerically that the f -mode is significantly perturbed in the presence of subsurface magnetic fields concentrated one to a few scale heights below the photosphere. The fanning (or increase in line width) of the f -mode was first reported in Singh et al. (2014) based on harmonically varying magnetic fields, but the associated mode strengthening was not emphasised there. Here we extended this by also investigating the effects of localised bipolar magnetic fields, resembling more realistically the active region precursors. The fanning effect of the f -mode is found to be weaker in the case of localised magnetic field concentrations. Motivated by the observational findings of Singh et al. (2016), where the strengthening of the f -mode at large wavenumbers prior to the emergence of about half a dozen active regions, we focus here primarily on the phenomenon of f -mode strengthening. In our numerical investigations reported here, the f -mode strength is found to be clearly larger at high horizontal wavenumbers in a variety of magnetic background states, compared to their corresponding hydrodynamic cases. While the fanning of the f -mode at large k x is more sensitive to the strength of the magnetic field at a given depth, the mode strength itself shows a dependence on the location of the magnetic concentration beneath the surface. Thus, it is indeed remarkable that the properties of the f -mode trace so effectively both the location and the strength of subsurface magnetic fields which are not yet directly visible at the photosphere.
We note that damping of the f -mode due to its resonant coupling with atmospheric magnetic fields or due to absorption caused by the sunspots is known (Cally and Bogdan 1997;Erdélyi 2006;Pintér et al. 2007), and also seen in observational study of Singh et al. (2016) after the emergence of active regions. But, the strengthening of the f -mode due to magnetic fields confined below the photosphere is a qualitatively different finding, which is numerically confirmed in the present work.
Let us now estimate the instrumental requirements to robustly probe and detect solar f -mode perturbations of the kind explored in this paper. In order to detect both fanning and mode-strength gain, according to figures 3 and 5,k x 2 needs to be reached. With γ H d ≈ 0.5 Mm and R = 700 Mm being the solar radius, this corresponds to spherical harmonic degrees 2800. However, indications of f -mode strengthening prior to active region formation could indeed be available even at somewhat smaller degrees, as was seen in the observational work of Singh et al. (2016) based on data from the HMI instrument, whose sensitivity drops beyond ≈ 2000. The precursor signal is expected to improve at larger wavenumbers. Therefore higher resolution data from other facilities and upcoming missions will be critical not only for establishing the connection between the solar f -mode and the subsurface magnetic fields, but also for providing deeper insights into the sunspot formation mechanism, which is still an outstanding topic in solar physics. More work is needed to understand the physical mechanism responsible for the strengthening of the f -mode and this will be attempted elsewhere. In particular, it would be useful to extend these models to explore effects of 3-dimensional magnetic configurations; see the review of Erdélyi (2006) for possible effects on mode damping and mode conversion.

Funding
Financial support from the Max Planck Princeton Centre for Plasma Physics (NS), Deutsche Forschungsgemeinschaft Heisenberg programme (grant No. KA 4825/1-1; PJK), the Academy of Finland ReSoLVE Center of Excellence (grant No. 307411; MJK, MR, PJK), the National Science Foundation Astrophysics and Astronomy Grant Program grant AST1615100, and the University of Colorado through its support of the George Ellery Hale visiting faculty appointment, is acknowledged.