BPS Boojums in N=2 supersymmetric gauge theories I

We study 1/4 Bogomol'nyi-Prasad-Sommerfield (BPS) composite solitons of vortex strings, domain walls and boojums in ${\cal N}=2$ supersymmetric Abelian gauge theories in four dimensions. We obtain solutions to the 1/4 BPS equations with the finite gauge coupling constant. To obtain numerical solutions for generic coupling constants, we construct globally correct approximate functions which allow us to easily find fixed points of gradient flow equations. We analytically/numerically confirm that the negative mass of a single boojum appearingat the endpoint of the vortex string on the logarithmically bent domain wall is equal to the half-mass of the 't Hooft-Polyakov monopole. We examine various configurations and clarify how the shape of the boojum depends on the coupling constants and moduli parameters. We also find analytic solutions to the 1/4 BPS equations for specific values of the coupling constants.


Introduction and summary
Topological solitons have long been appreciated in various fields in modern physics such as string theory, field theory, cosmology, nuclear physics and condensed matter physics. They often appear with spontaneously broken symmetries which support the stability of the solutions with conserved topological charges. Typical examples are Nielsen-Olesen vortex string in the Abelian-Higgs model [1] and 't Hooft-Polyakov magnetic monopole in SU (2) Yang-Mills theory [2,3]. The topological charges are associated with the homotopy groups π 1 (U (1)) and π 2 (SU (2)/U (1)), respectively. Domain walls also often appear in models with degenerate and discrete vacua which are characterized by the homotopy group π 0 (M vac ), where M vac is a vacuum manifold.
In addition to these elementary solitons, there are also composite solitons, which include two or more different kinds of elementary solitons. They are sometimes called D-brane solitons or field theoretical D-branes, because they share many properties with D-branes in superstring theories [4]. Among various configurations, in this paper, we are especially interested in composite solitons, where vortex strings attach to domain walls. Such configurations are natural field theory counterparts to the Dbranes on which fundamental strings end in superstring theory. These configurations have been studied for a long time. A simple field theoretical model possessing the composite soliton was provided in Ref. [5]. Numerical analysis based on the specific scalar field theory was performed in [6]. The existence of such configurations in N = 1 supersymmetric (SUSY) QCD was found in [7], and inspired by this result, qualitative explanation of existence of a vortex-string ending on a domain wall was shown in [8,9]. It was found that annihilation of domain wall and anti-domain wall produces stable solitons of lower dimensionality [10] in a modified model of N = 1 SUSY QCD [7].
Much progress has been done within SUSY gauge theories in the past decade. In N = 2 SQED, seen as a low energy effective theory of N = 2 SQCD with SU (2) gauge group with N F = 2 hypermultiplets, perturbed by a small mass term for the adjoint scalar field, the 1/4 BPS equations for the Abelian vortex strings ending on the domain walls were first derived in [11]. Endpoints of the vortex strings on the domain wall were identified with electric particles in a low energy effective theory of the domain wall. Then the model was extended to N = 2 SU (2) × U (1) SQCD with N F = 4 hypermultiplets [12] and the 1/4 BPS equations including the non-Abelian vortex strings [13,14] were found. A low energy effective theory on the composite domain walls was studied and it was found that the endpoints of the non-Abelian vortex strings play the role of a non-Abelian charge coupled with localized dual non-Abelian gauge fields at linear level of small fluctuations [12], see also related works [15,16,17,18,19].
Almost in the same period, the different kind of the 1/4 BPS composite con-figurations of the non-Abelian vortex and the magnetic monopoles were found in N = 2 U (N C ) SQCD with N F = N C flavors [14,20,21,22,23], see also recent related works [24,25,26,27,28]. Afterward, all these kinds of topological solitons, the vortex strings, domain walls, and the magnetic monopoles were found to coexist as the 1/4 BPS composite states in N = 2 U (N C ) SQCD with N F ≥ N C flavors. The most generic 1/4 BPS equations for these were found in [29]. The BPS equations are a set of first order differential equations for the gauge fields A a µ , the adjoint scalar fields Σ a (a = 1, 2, · · · , N 2 C ), and N C × N F squark fields H andĤ, which depends on three spatial coordinates x k (k = 1, 2, 3). Since, except for several simplest configurations, the generic solutions have no spatial symmetries, it is quite hard to solve the BPS equations, even though they are first order differential equations. To make things worse, no analytic solutions to the 1/4 BPS equations have been found in the model with finite gauge coupling constants. In order to overcome these difficulties a powerful technique, so-called moduli matrix formalism, was invented for the 1/2 BPS domain walls [30,31], 1/2 BPS vortex strings [32] and for the full 1/4 BPS equations [29]. It is so powerful that not only the dimension of the moduli space but also all the information on the moduli parameters for all topological sectors are easily exhausted. Notably, the moduli matrix formalism provides all exact solutions in the strong gauge coupling limit where the model reduces to the massive nonlinear sigma model whose target space is a cotangent bundle over Grassmannian manifold Gr N F ,N C SU (N F )/SU (N C ) × SU (N F − N C ) × U (1) [29]. There is another advantage of the moduli matrix formalism, though less emphasized in the literature, that the complicated first order differential equations for 2N 2 C fields in the vector multiplet and 2N C N F fields in the hypermultiplet (Ĥ = 0) are converted into N C (N C + 1)/2 second order equations, collectively called the master equation for real positive N C × N C matrix field Ω [29], where g is the U (N C ) gauge coupling constant, v 2 is the FI parameter and Ω 0 stands for a source term, which can be chosen quite freely up to certain rules [29]. Even in the case N F = N C , for which the moduli matrix formalism gives a minimal benefit, the degrees of freedom reduces by half. Details of the moduli matrix formalism are summarized in the review paper [33]. Also, there are good reviews covering many studies of the topological solitons in N = 2 SQCD [34,35,36,37]. While many properties are shared by BPS topological solitons in field theory and D-branes in string theory, there is a BPS object, which is inherently field-theoretical and has no analog in string theory. It appears at a junction point where two different kinds of topological solitons meet. Interestingly, it has a negative BPS mass and can be interpreted as a binding energy. The first example of this kind of object was found as the domain wall junction [38,39,40,41,42,43]. When a vortex string ends on a domain wall, a different type of BPS object with negative mass appears [29]. Like the domain wall junctions, it was shown that the negative BPS mass is nothing but the binding energy of the vortex string and the domain wall [44]. In particle physics context it is now called the boojum [44] because a similar configuration, coined Boojum 1 by Mermin, is known in the context of 3 He superfluid [45,46]. Recently, the boojums are getting more popular in various fields. For example, it has been studied in the 2 component Bose-Einstein condensates [47,48], and also in the dense QCD [49,50]. The negative BPS binding energy was also found for intersections of vortex sheets in 5 dimensions [21,51,52]. Typically, these negative BPS masses appear in Abelian-Higgs model and these three different binding energies in different dimensions can be reasonably understood through a descent relations given by Kaluza-Klein dimensional reductions [53].
Our main aim in this paper is to study the 1/4 BPS boojums in details in N = 2 SUSY QED with N F ≥ 2 flavors in the presence of the Fayet-Iliopoulos term. If one desires, one can think of this theory as a low energy effective theory of the mass perturbed SU (2) SQCD as considered in [11]. As mentioned above, fundamental results such as the derivation of the BPS equations, the BPS masses and the structure of the moduli space were already done [29,33], and the interaction rules of various solitons and some qualitative features were clarified [44]. However, we would like to point out that these understandings still remain at the qualitative level. This is because neither analytic nor numerical solutions to the 1/4 BPS equations have been obtained. Although all exact solutions to the 1/4 BPS equations were obtained [29], it was done only in the strong gauge coupling limit where the BPS boojum masses becomes zero and the lump strings in the bulk are singular. In order to fill this hole, we will provide both numerical solutions for generic case and analytical solutions for particular cases with the gauge coupling constant kept finite. To this end, we are greatly helped by the moduli matrix formalism [29,33], which reduces the complicated 1/4 BPS equations to a mere second order equation for a single real scalar field u(x k ) 1 2g 2 ∂ 2 k u = v 2 1 − Ω 0 e −u , (k = 1, 2, 3), (1.2) which is the Abelian version of the master equation given in Eq. (1.1) with identification Ω = e u . In this paper, we will, in particular explain how to solve this Abelian master equation. At this point, it is worth mentioning that the two-dimensional version of Eq. (1.2) is the well-known Taubes equation for the BPS vortices in the Abelian-Higgs theory [54]. Unfortunately, no analytic solutions are known for the Taubes equation in R 2 . This is in contrast to the familiar BPS Yang-Mills instantons, for which the well-known Atiyah-Drinfeld-Hitchin-Manin construction is established [55]. However, Taubes equation on the hyperbolic plane H 2 of curvature −1/2 is slightly modified from Eq. (1.2), such that the first term in the parentheses on the right-hand side vanishes and the equation becomes integrable [56]. This is a consequence of the fact that these vortices are obtained as a dimensional reduction from SO(3)-symmetric Yang-Mills instantons on R 4 [56]. Recent progress for integrable vortex models are found in [57,58,59,60,61,62,63]. In short, no analytic solutions have been found for the Taubes equation in R 2 , much less for the master equation (1.2) in R 3 . Thus, it seems almost hopeless to find analytic solutions to the master equations (1.2). Surprisingly, we will provide some exact solutions in both R 2 and R 3 .
A technical but important result of this paper is finding a global approximation to solutions of the master equation (1.2), which supports almost all other results in this paper. In general, solving the Cauchy boundary problem in R 3 for configuration of identical topological solitons, such as parallel domain walls or parallel vortex strings, is not difficult because finding appropriate boundary conditions offers no complications. However, when two or more topological solitons of a different kind coexist, giving an appropriate boundary condition is not always a straightforward task. A prototype example is a vortex string attached to a domain wall from one side. The vortex string has codimension two while the domain wall has codimension one. As was mentioned above, the domain wall is pulled and it logarithmically bends toward spatial infinity along the string axis. Thus, in order to solve the Cauchy boundary problem for this kind of configuration, we have to have an asymptotic solution for bent domain wall. We will provide a simple but very generic method of constructing globally correct approximate solutions from the solutions of constituent isolated solitons. Such global approximations are very close to the exact solution almost everywhere except around vicinity of the junction points. Furthermore, they are regular everywhere unlike the solutions in the strong gauge coupling limit, where vortex-strings develop singular cores. We will use these global approximations not only for determining the boundary condition of the master equation (1.2) but also for a suitable initial function consistent with the boundary condition to solve Eq. (1.2) by the gradient flow (imaginary time relaxation) method. The closer initial function is to the true solution, the faster the gradient flow converges. Therefore, the global approximate solution is useful for numerical works. Whatever the solutions for the constituent solitons are, numerical or analytical, our method works very efficiently. Namely, with the global approximations, we can solve the Cauchy problem for any kind of configuration in the theory with arbitrary coupling constants.
As a direct consequence of solving the BPS equations for the finite gauge coupling, we reveal shape of the boojums, see Fig. 4. So far, only a schematic picture such as a simple hemisphere for the boojum have been given in the literature. We will also show how the shape of the boojums is modified when the coupling constants of the model are changed. It is also interesting to observe how they deform when multiple boojums coalesce. Do they behave similarly to the BPS magnetic monopoles, that it that two separate balls collide and form a donut [64]? Our numerical computations show that the boojums behave very differently from the monopoles, see Figs. 12, 13 and 21. Furthermore, the global approximation developed in this paper settles down a problem raised in [65]: Based on approximate solutions it was pointed out that there is an ambiguity in the definition of the Boojum mass, which stems from the ambiguity of the definition of the geometric parameters such as domain wall area and vortex string length caused by logarithmic bending of the domain wall. As mentioned above, the boojum mass was originally computed in [44] and it was shown that it is a negative half mass of 't Hooft-Polyakov monopole without any ambiguity. To reach this value, they simplified their calculation by taking a mean value of the scalar field in the vector multiplet at a cross section of logarithmically bend domain wall [44]. Although the result -the negative half of 't Hooft-Polyakov monopole -seems to be plausible, the validity for taking the average is not very clear because the cross section is exponentially large far away from the junction point. We will recalculate the boojum mass using the global approximation and confirm that the calculation in [44] is correct without any ambiguity. In [65], they also considered the case of two vortex-strings attached from the both sides of the domain wall, ending at different points (see Fig. 14). According to [65], this setting raises another problem about localization of the binding energy. Two possibilities are considered: One is that the binding energy locates around the junction points only and the other is that the binding energy is also localized half way between the strings endings. Our numerical computation confirms the former scenario, as is clearly visible in Fig. 10.
We have another progress. We have stressed that no analytic solutions have been found for the master equation (1.2) so far. The situation is the same even for a much simpler case of the Taubes equation, Eq. (1.2) in R 2 , for the local vortices. Thus, one might think that there is no hope to find analytic solutions to the master equations (1.2) either in R 2 or R 3 . Contrary to expectations, we will obtain several analytic exact solutions to the Taubes equation for specific semi-local vortices in R 2 . Furthermore, combining this with the known exact solutions to the master equation in R 1 for the domain walls [66] in a similar way with which we made the global approximate solution, we will, surprisingly, succeed in constructing analytic solutions to the master equation (1.2) in R 3 . In connection with these results, in Appendix B we also develop accurate approximation to the individual vortex string and domain wall solutions, which are useful to quickly obtain the global approximation for the 1/4 BPS solutions.
This paper is only a first half of a two-paper miniseries where we present our results on 1/4 BPS solitons. Here, we primarily deal with the technical issues connected with solving 1/4 BPS equations and we investigate elementary properties of the Boojum, such as its mass and shape. In the second paper [67], we collect our results about 1/4 BPS solitons which are more in-depth. First, by observing the two-dimensional spreading of magnetic field lines from the boojum inside a thick domain wall and matching it to the Coulomb law for a point magnetic charge at an asymptotic distance, we cement the notion of a Boojum being a confined fractional magnetic monopole. We further pursue this analogy by defining a "magnetic scalar potential", which we identify as a solution to the vortex part of the master equation. A novel solution of a semi-local boojum, which arise when a semi-local string with a size moduli attaches to a domain wall, is also investigated. In addition, we study configurations with an infinite number of vortex-strings aligned on one or both sides of the domain wall and quantify the ability of such configurations to store magnetic charge as "magnetic capacitors". Dyonic extensions of 1/4 BPS solitons are also pursued. Most notably, we discover that the positive and negative electric charge density is stored on opposite skins of the domain wall, making a thick domain wall an "electric capacitor". Finally, we also study 1/4 BPS configuration from the viewpoint of the low energy effective action, the Nambu-Goto action and the DBI action, for the domain wall. This paper is organized as follows. In Section 2, we review our model, which is N = 2 supersymmetric U (1) gauge theory in (3 + 1)-dimensions coupled to N F hypermultiplets. We derive the BPS equations for the boojum configuration, preserving 1/4 supersymmetry and shows that they reduce to one differential equation (1.2). In Section 3, we deal with the axially symmetric configuration, namely one vortex string ending on the domain wall. We show all the basic properties of the boojum there. We also calculate the boojum mass there. Then, we demonstrate more complicated solutions, where multiple vortex strings end on a single domain wall in Section 4. We study similar configuration with two or more domain walls in Appendix A. We summarize all global approximations for 1/4 BPS solutions in section 5. These are complemented by accurate approximations to the ANO vortex and the domain wall which we develop in Appendix B. Several exact solutions to the 1/2 and 1/4 BPS equations in presented Section 6. A brief discussion of the future work is given in Section 7.

Abelian vortex-wall system
Let us consider N = 2 supersymmetric U (1) gauge theory in (3+1)-dimensions with 2N F complex scalar fields in the charged hypermultiplets. We assemble them into a row vector H ≡ (H 1 , H 2 , . . . , H N F ) andĤ † ≡ (Ĥ 1 * ,Ĥ 2 * , . . . ,Ĥ N F * ). The vector multiplet includes the photon A µ and a real scalar field σ. The bosonic Lagrangian is given as where g is a gauge coupling constant, M is a real diagonal matrix and v is the Fayet-Illiopoulos D-term. Without loss of generality we can take M to be traceless, namely N F A=1 m A = 0, 2 and align the masses as m A > m A+1 . SinceH will play no role, we will setH = 0 in the rest of this paper. Throughout this paper, we use the conventions: Since our model is supersymmetric, we can reduce the full equations of motion into BPS equations, which are first order differential equations. Solutions to the BPS equations are called BPS solitons. These field configurations possess many special properties, such as saturating the energy bound in a given topological sector and preserving a fraction of supercharges in the supersymmetric theory. These facts simplify construction and discussion of topological solutions significantly. Of course, one can consider more general values of coupling constants without qualitatively changing the results. We will focus on the BPS solutions because the generic case would involve much-complicated analysis for solving the full equations of motion.
In the absence of the mass matrix M , the Lagrangian (2.1) is invariant under SU (N F ) flavour transformation of Higgs fields H → HU , U ∈ SU (N F ). The nondegenerate masses in M explicitly break this down to U (1) N F −1 , which we from now on assume to be the case unless stated otherwise.
There are N F discrete vacua with representative field values We can act on each of these configurations by N F global U (1) transformations, which make each vacuum topologically a circle. Thus, the space of all vacuum configurations -the manifold V -is isomorphic to N F -torus V ∼ S 1 × S 1 × . . . × S 1 . A direct consequence of this fact is that the model (2.1) has a rich spectrum of topological excitations [20,29,44]. Non-trivial homotopy groups π 0 (V) = Z N F and π 1 (V) = Z × Z × . . . × Z give rise to a large variety of domains walls and vortex strings. These solitons are 1/2 BPS solitons, as they preserve half of the supercharges if the model is extended to have supersymmetry. In addition, there are also lots of composite solitons, where vortex strings attach to the domain walls. These configurations preserve 1/4 of supercharges and hence are known as 1/4 BPS solitons, which we are interested in.

1/4 BPS state
Let us now construct the 1/4 BPS solitons. We will arrange vortex strings to be parallel to the x 3 axis and the domain walls to be perpendicular to the x 3 axis. We let A 0 = 0 and consider static configurations, so that the Higgs field H, the scalar σ and the gauge field A k (k = 1, 2, 3) are independent on x 0 . The Bogomol'nyi bound is then obtained as 3 where ε klm is a totally antisymmetric symbol in three dimensions (ε 123 = 1) and j k are non-topological currents defined by The energy density (2.11) saturates the so-called Bogomol'nyi bound with The following useful identity holds: if the following 1/4 BPS equations are satisfied. It is easy to verify that these equations are compatible with the equations of motion (2.7)-(2.9) for all values of parameters η 2 = ξ 2 = 1, where ξ = (−1)1 labels (anti-)vortices and η = (−1)1 denotes (anti-)walls. T W and T S , the domain wall and the vortex string energy density respectively, are positive definite. T B is the so-called boojum energy density, which is interpreted as binding energy of vortex string attached to the domain wall, since it is negative irrespective the signs of η and ξ [44,65]. The total energy of 1/4 BPS soliton is obtained upon the space integration and it consists of three parts Summing up all the elementary domain walls and vortex strings, we have where we have denoted ∆m = [σ] x 3 =+∞ x 3 =−∞ and k ∈ Z stands for the number of vortex strings. T B has been also calculated in [29,44,65]. But, there is a discussion about determining the mass of a single boojum which appears as a junction of the semiinfinite string ending on a logarithmically bent domain wall. In what follows, we will confirm without doubt that the boojum mass even for bent domain walls is given by the following formula [44] where the sum is taken for all the junctions of domain walls and vortex strings in the solution under consideration.

The moduli matrix formalism
With the use of the moduli matrix approach [30,29,33] where H 0 (z) is the so-called the moduli matrix which is holomorphic in a complex coordinate z ≡ x 1 + iξx 2 . Our notation is ∂z = (∂ 1 + iξ∂ 2 )/2. By using the U (1) gauge transformation, we fix u = u(z,z, x 3 ) to be real. Then we have A 3 = 0. The last equation (2.19) turns into the master equation Now, all fields can be expressed in terms of u as follows The energy densities are also written as The non-topological current j k given in Eqs. (2.12) and (2.13) can be rewritten in the following expression by using the BPS equations Thus, we also have Collecting all pieces, the total energy density is given by Thus, the scalar function u determines everything. The moduli matrix formalism is an easy tool for describing the 1/4 BPS solitons. Indeed, one can construct any kind of configurations by preparing appropriately the moduli matrix H 0 (z) = (h 1 (z), h 2 (z), · · · , h N F (z)). As a simple example, let us consider the case that only h 1 and h 2 are non-zero. The vacuum i is given by h j = δ ij (j = 1, 2). If they are both non zero constants, the corresponding configuration has the domain wall separating 1 and 2 vacua. Since we assume m A > m A+1 , the vacuum 1 appears on the left hand side of the domain wall while the vacuum 2 on the right hand side. When one wants to put N 1 vortex strings in the vacua 1 , one only needs to prepare h 1 = P 1 (z) where P 1 (z) stands for a polynomial with degree N 1 . This can be straightforwardly extended for generic H 0 (z).
Hence, the moduli matrix formalism allows us to handle the complicated 1/4 BPS solitons. One can even deal with the solutions belonging to different topological sectors simultaneously. It is nice that all information about the moduli parameters are included in H 0 . For these reasons, the moduli matrix formalism has been frequently used in the literature since its inception. Nevertheless, a finishing touch has been omitted in the sense that no analytic or numerical solutions have been constructed for the master equation (2.25) except for the strong gauge coupling limit g 2 → ∞ [30] as was mentioned in the Introduction. In the literature, only schematic pictures in the weak/strong gauge coupling region has been given.This lack of solutions also caused some confusion about the definition of the boojum [65], which we will mention later.

Axially symmetric solutions
In this section, we explain how we numerically solve the master equation (2.25). Throughout this section, we will concentrate on N F = 2 case but everything addressed here can be straightforwardly extended to generic case. For ease of notation, we will use the following dimensionless coordinates and mass The dimensionless fields are similarly defined bỹ We will also use the dimensionless magnetic fields Then, the dimensionless energy densityẼ is defined by The relations to the original values are given as In what follows, we will not distinguish x k andx k , unless stated otherwise. An exception is the mass: we will use the notationm in order not to forget that we are using the dimensionless variables. Then the master equation (2.25) is expressed as If the moduli matrix is constant, we have flat domain walls and the function u = u W (x 3 ) depend only on x 3 coordinate. The reduced master equation reads Similarly, when the moduli matrix has only one non-zero component, say H 0 = (P n (z), 0, · · · , 0), the configuration has n vortex strings in the vacuum 1 . One can eliminate the x 3 dependence by defining Then the master equation reduces to

Gradient flow
Now we are ready to solve the master equation Eq. (3.14) numerically. But, instead of solving it directly, we will solve the so-called gradient flow equation for U = U (x k , t), using appropriate initial configuration U (x k , 0) and boundary conditions. Normally, U (x i , t) rapidly converges to a static function. When ∂ t U (x k , t) becomes negligibly small at a sufficiently large t, we regard U (x k , t) as a solution u(x k ) of the original master equation, The final state obtained in this way depends on the initial configuration and the boundary conditions. Thus, the necessary step to successfully solve Eq. (3.18) is to provide an initial configuration compatible with the boundary conditions.

Domain wall
Solving the gradient flow equation is not difficult especially for a single type of solitons like the flat domain walls or the straight vortex strings. For the domain walls perpendicular to the x 3 axis, the master equation reduces to Eq. (3.15). The solution u W (x 3 ) can be obtained by solving the corresponding gradient flow equation A convenient choice of the initial configuration and the boundary conditions at The true solution can be obtained as the asymptotic function u W ( Several numerical solutions are shown in Fig. 1. There are two qualitatively different domain walls according tom [11]. In the strong gauge coupling limit, m 1, the domain wall has a simple structure as shown in the left panel of Fig. 1. The traverse size of the domain wall isd W = 2/m. On the other hand, in the weak coupling region,m 1, the domain walls consists of two layers [11]: the fundamental Higgs field decays in the outer thin layer whose transverse size is of order 1 (in the dimensionless unit), whereas the singlet scalar fieldσ interpolatesm/2 and −m/2 through the inner fat layer whose transverse size is aboutd W = 2m. In the latter case, the inside of the domain wall remains almost the Coulomb phase since the charged fundamental fields are exponentially small there.

Vortex string
The vortex strings can be obtained similarly. As an example, let us consider n vortex strings in the 1 vacuum. It is given by the moduli matrix H 0 = (P n (z), 0) with P n (z) being the n-th order polynomial. This yields Ω 0 = |P n | 2 e ηmx 3 . Then we define This yieldsσ = ηm. The master equation to U S for the vortex strings reduces to for the initial and boundary conditions is where, of course, we assume |L 1,2 | 1 so that e −ρ 2 can be neglected. The exponential term appearing in Eq. (3.25) is there in order to avoid singular behavior at the center of the vortex strings in the initial configuration. It is exponentially small at the boundary x 1 = ±L 1 and x 2 = ±L 2 , so that the initial configuration does not contradict to the boundary conditions. As before, the true vortex solution u S for Eq. (3.17) can be obtained as the asymptotic function u S (x 1 , The single vortex string (P 1 = z) leads to an axially symmetric configuration U S (ρ, t) satisfying the master equation This should be solved with the boundary conditions The asymptotic behavior of u S (ρ) = lim t→∞ U S (ρ, t) is The numerical solution for a single vortex string is given in Fig. 2.
x 3 x 1 Figure 3: The vortex trumpet: vortex string ending on a logarithmically curved domain wall. The asymptotic behavior of the function u(ρ, x 3 ) is shown on each part of the boundary (the 'long cylinder').

Vortex-string attached to a domain wall
Let us now address how to get the simplest 1/4 BPS solutions defined by the moduli matrix H 0 = (z, 1). The corresponding gradient flow equation is given by A significant difference from the homogeneous solitons is that the boundaries of the 3D box parallel to the vortex string intersect with the curved domain wall and one of the boundary parallel to the domain wall intersects with the vortex string, see Fig. 3. Thus, we have to figure out an appropriate boundary condition consistent with the curved domain wall. Let us first look at the region far from the vortex string axis, namely ρ 1. Since u depends on x 1 and x 2 only through log ρ, we can drop ∂ 2 ρ + (1/ρ)∂ ρ term from the master equation This can be rewritten as This is nothing but the master equation for the domain wall if we regard u − log ρ The shift (log ρ)/m in the argument of u W reflects the fact that the domain wall asymptotically bends as Next, going far from the domain wall, u approaches to the vortex string u → mx 3 + u S for x 3 0 while becomes close to the vacuum u → −mx 3 for x 3 0.
A suitable function possessing these asymptotic behaviors can be obtained by replacing ρ by e u S /2 in the asymptotic function u given in Eq. (3.33). This is our choice for the initial configuration to solve the gradient flow , which is desired asymptotic behavior of Eq. (3.35). The boundary conditions consistent with the initial configuration are Note that the first condition is ensured by ∂ ρ u S (ρ = 0) = 0. The initial configuration (3.36) gives us an important information about the domain wall's position for whole ρ as which is consistent with the asymptotic behavior given in Eq. (3.34). In the strong gauge coupling limit, we have the exact solution u S = log ρ 2 , so the asymptotic relation (3.34) becomes exact for whole ρ. As shown in the left panel of Fig. 2, the domain wall is smooth everywhere in the finite gauge coupling, but it gets logarithmically sharp at the origin in the strong gauge coupling limit.

Boojum at different gauge couplings
Having the appropriate initial configuration U, we now solve the gradient flow equation (3.30). We solve it by using a standard finite difference method with the Crank-Nicolson scheme. A typical solution withm = 1 is shown in Fig. 4. One can clearly see that the straight vortex string with a finite diameter ends on the logarithmically curved domain wall. This is an obvious contrast to previously obtained solutions with singular vortex strings [68,29]. Furthermore, most importantly, our solution has the boojum which contributes to the total energy (the boojum disappears in the strong coupling limit since the boojum charge is proportional to 1/g 2 ). With our numerical solution at hand, we figure out the correct shape of the boojum as shown in Fig. 4. So far, the shape of boojum has been only schematically realized. For example, in Ref. [44], the boojum was expressed as just a hemisphere which is little too simple compared to our numerical solution. The third advantage to have the numerical solution for a finite gauge coupling is that we can clearly see the distribution of the magnetic flux which enters from the vortex string into the domain wall. In Fig. 4 we show several lines of the magnetic field (blue solid curves). They form a flux tube inside the vortex string and radially spread out once they enter the domain wall. This picture can correctly be obtained only for the finite gauge coupling since electromagnetic interaction disappears in the strong coupling limit. The numerical solutions also shed light on both qualitative and quantitative difference between the 1/4 BPS configurations in the strong (m 1) and the weak (m 1) gauge coupling regions. Especially, the vicinity of a junction point of the domain wall and vortex string is hard to correctly understand without numerical methods. We show the several solutions in details; the profiles of the fields, the topological charge densities, and the energy densities are displayed in Since we have rescaled out gv dependence, the asymptotic vortex string configuration is common for all the cases. This can be seen in regions near the top-left edges of the middle panels of Figs. 5 -7 in which the topological charge densities of the vortex stringT S are shown. The bottom-middle panels show the magnitude of the magnetic fields which are slightly different from the string topological charge densityT S .
The three panels in the first row of Figs. 5-7 showH 1 ,H 2 andσ. The fieldH 1 develops the non-zero VEV in the upper region (x 3 0). Namely, the light yellow regions of the top-left panels correspond to the vacuum 1 . Similarly, the yellow parts of the top-middle panels show the region of the vacuum 2 . Since the phase ofH 1 has the winding number 1 for the vortex string,H 1 vanishes at the core of the vortex string (on the top-left edge (ρ = 0, x 3 0)). Since neitherH 1 norH 2 develops non-zero VEV there, the U (1) gauge symmetry is recovered.
Shape of the bent domain wall can be seen in the middle-left panels where the topological charge densityT W is plotted. The yellow-line region ofT W is consistent with the yellow-blue transit region ofσ in the top-right panels. The domain wall strongly bends for the strong coupling (m 1) while it seems to be almost flat for the weak coupling (m 1). This difference comes from the wall to string tension ratio, which can be found by comparing the scales of the color bars of the middle-left and middle-middle panels.
The domain wall tends to strongly bend when the tension of the vortex string is larger than that of the domain wall. This can be also understood from Eq. (3.34), which tells that the asymptotic wall curve is x 3 − log ρ 1/m . For example, x 3 (ρ = 10) = − log 10 for the case ofm = 1. To get the same amount of bending in the case ofm = 10, one needs to go to ρ = 10 10 ! This is the reason why the domain wall seems almost flat form = 10 in the middle-left panel of Fig. 7. It is also seen that the domain wall at ρ = 0 is not singular for anym. In contrast to the inside of the vortex string, there are no topological reasons for the domain wall core to be in the unbroken phase. Indeed, bothH 1 andH 2 are non-zero inside the domain wall with m = 1/5 and 1. This is consistent with the single domain wall solutions given in Fig. 1. Only when the gauge coupling is very weak (m 1), the domain wall core develops the fat inner layer of the unbroken phase. While the diameter of the vortex string is fixed to be of order one, width of the domain wall varies asm is changed, see Fig. 1. The domain wall is sharpest whenm 1, while it becomes much fatter than the vortex string for both weak and strong gauge coupling limits, see the bottom-left panels showing the full energy density including the surface term.
Our prime interest is to understand the junction point of the domain wall and the vortex string. Especially in the weak coupling limit, since the cores of the vortex string and domain wall are in the unbroken phase. The top-left and top-middle panels of Fig. 7 clearly show that the two unbroken regions are just smoothly connected. We are also interested in the boojum, the binding energy appearing at the junction point. The middle-right panels exhibit a variety of shapes of boojum depending oñ m. It is basically a smoothed cone which becomes thinner form 1 and fatter for m 1. Finally, in the bottom-right panels we plot the mass square of the gauge field defined by The gauge symmetry is recovered in the deep-blue region, while is strongly broken in the light-yellow region. In addition, we also plot flow of the magnetic field (magnetic force lines). In general, the flow tends to be concentrated into the region where the gauge symmetry is unbroken or weakly broken. In the case of the strong gauge coupling (m = 1/5 given in Fig. 5), the incoming magnetic force lines through the vortex string are forced to bend when they encounter the domain wall. They scatter and then spread out along the domain wall because the breaking of the gauge symmetry is milder inside the domain wall. However, the lines running near the core of the vortex string are pushed out of the domain wall. This happens because the energy cost to bend the lines all the way into the core of domain wall is larger than that to have mildly bending lines going through the bulk. In contrast, in the intermediate case withm = 1 given in Fig. 6, the magnetic force lines go along the domain wall. This is because the mass of the gauge field inside the domain wall is sufficiently small so that the magnetic force lines bend and tend to go inside the domain wall. The weak coupling case withm = 10 shown in Fig. 7

Boojum mass
Finally, let us compute the tensions of the topological objects: The domain wall tension is everywhere constant in the x 1 -x 2 plane and it is given as The string tension is a constant along the x 3 directioñ where we used the asymptotic behavior given in Eq. (3.33) and u W (x 3 ) →mx 3 for x 3 0. In terms of the original variables, these are rewritten as One should be careful that the (dimensionless) string tension is always 4π even in the vacuum 2 side of the domain wall (x 3 0) where we do not put the vortex string. This is simply a consequence of the flux conservation. The magnetic flux cannot invade the vacuum 2 , so that it flows along the surface of the domain wall which logarithmically bends. Therefore, regardless of x 3 , we always have a total flux of 4π. In other words, the vortex string is not chopped by the domain wall but exponentially inflates. This observation can also be justified by seeing the profile of the scalar field H 1 . The amplitude of the scalar field whose phase has the winding number must necessarily vanish somewhere by a topological reason. Therefore, we can think of the deep-blue region in the top-left panels of Figs. 5 -7 as the inside the vortex trumpet.
The total magnetic flux outgoing through a sufficiently long cylinder (the right side is at x 3 R m and the left side is at x 3 L −m) of sufficiently large radius R ( 1) surrounding the vortex trumpet, see Fig. 3, can be calculated as follows where we used Eqs. (3.22) and (3.33). The incoming magnetic flux passing the right side of the cylinder can be similarly obtained as where we have used the boundary conditions, ∂ ρ u S → 0 for ρ → 0 and u S → log ρ 2 for ρ → ∞. As expected, we haveΦ C +Φ R = 0 andT S = −Φ R due to the flux conservation. This conservation implies there is no magnetic flux passing the left side of the cylinder, namelyΦ L = 0. In order to express this in a conventional way, we need to use the original variables and rescale A µ → gA µ . Then we have Similarly, the boojum mass is calculated as follows.
where we usedσ →m/2 for x 3 = x 3 R m and the following identity Note that the above result reflects an accidental Z 2 symmetry for our specific choice of the mass matrixM = diag(m/2, −m/2). For generic case withM = diag(m 1 ,m 2 ) (m 1 >m 2 ), the asymptotic behavior of u at a larger ρ becomes u → u W ( . It is straightforward to verify the following identity cylinder dS i ijk σF jk = 4π (m 1 +m 2 ) . (3.49) Thus, we find the generic form In terms of the original variables, this can be expressed as This is consistent with the previously obtained value in Ref. [44].
As a closing comment for this subsection, we would like to emphasize that we were able to rigorously calculate the charges thanks to the asymptotic behavior given in Eq. (3.33). In particular, we calculated the charges by surrounding them with the finite size cylinder as shown in Fig. 3. Since we have not touched the spatial infinity, the above computations do not suffer from any complications coming from logarithmic bending at all. In the literature, several attempts to avoid this problem have been done. In Ref. [29], the mass of two coaxial boojums was calculated for the flat domain wall on which one coaxial vortex string ends on from both sides and found 2T B = −4π(m 1 − m 2 )/g 2 . Since the domain wall is flat in this case, there are no anxieties. Then, in the proceeding work [65], it was proposed that the mass of the single boojum can be obtained by just dividing the mass of the two coaxial boojums by 2. In the same paper [65], they also discussed that there is an uncertainty for determining the mass of the single boojum in the case of logarithmically bent domain wall due to a geometrical ambiguity. Ref. [44] is the first work which showed T B = −2π(m 1 − m 2 )/g 2 for the isolated single boojum in the logarithmically curved domain wall. They reached this value by integrating the mass density T B at under the assumption that no flux goes through the boundary transverse to the vortex string axis. The first term unambiguously equals to −2πm 1 /g 2 . The second term is somewhat tricky since one needs to consider the intersection of the logarithmically bent domain wall with the boundary at x 3 = −∞. In Ref. [44], σ was simply replaced by a mean value (m 1 + m 2 )/2 and made use of the flux conservation, then the second term was found to be −π(m 1 +m 2 )/g 2 . Summing up the two contribution gives us T B = −2π(m 1 −m 2 )/g 2 . This way of computing the boojum mass is plausible, but it might be better to confirm it by a more rigorous way. We believe that our calculations above contribute to this point. For further confirmation, we perform the integrations numerically and find −T B /8πm = 0.976, 0.996, 1.010, 0.997 form = 1, 10, 20, 30, respectively. In conclusion, we agree with [44] that the boojum mass can be cleanly separated from the semi-infinite vortex string and the logarithmically bent domain wall.

Non-axially symmetric solutions
In the previous section, we have concentrated on studying the 1/4 BPS states which are all axially symmetric about the vortex-string axis. This symmetry reduces the problem from being three-dimensional to two-dimensional. In this section we will show many numerical solutions, which require full three-dimensional treatment. We can again use the moduli matrix formalism introduced in Sec. Let us first look at the 1/4 BPS solutions in N F = 2 case withM = (m/2, −m/2), which has two distinct vacua 1 and 2 , and one domain wall separating those vacua. The moduli matrix for n 1 (n 2 ) vortex strings in the vacuum 1 ( 2 ) is given by where P n (z) stands for an arbitrary polynomial of n-th power in z. For this moduli matrix, we solve the master equation (3.14) written in terms of the dimensionless parameters (3.1). As before, we upgrade the master equation to the gradient flow (3.18) for the real scalar function U (x k , t), We need an appropriate initial function for solving this. It is given by where u W (x 3 ) is the domain wall solution to the master equation (3.15) and u (n) S (x 1 , x 2 ) is the n vortex string solution to the master equation (3.17). One can easily check that this correctly reproduces the initial functions given in Sec. 3 for n 1,2 = 0 or 1.
Two vortex strings in the first vacuum 1 ending on the domain wall can be generated by The vortex strings are asymptotically parallel and located at z = ±L. In Fig. 8 we show a numerical solution with L = 6 form = 1. The gray surfaces in the panels (a1) and (a2) on Fig. 8 are the energy density isosurfaces on whichT W +T S +T B takes one half of its maximum. Figure 8: The plots show the energy density isosurfaces of two vortices ending on one wall (a1, a2), where the blue and the red curves show magnetic flux, the wall energy density (b1), the vortex energy density (b2), the boojum energy density (b3) and the total energy density (b4) with the distance between two vortices L = 6. The mass is set to bem = 1. The red lumps in the panels (a1) and (a2) of Fig. 8 show the boojum density isosurface on whichT B takes one-half of its minimum (a negative value) and the red and blue curves correspond to the magnetic force lines going through the vortex strings.
In the panels (b1), (b2), and (b3) we show the topological charge densitiesT W ,T S , andT B on the cross section which passes the centers of vortex strings, respectively. The panel (b4) of Fig. 8 depicts the total energy densityT including the surface terms on the same cross section. One can easily change the distance of two vortex strings by varying L in Eq. (4.4). The numerical solutions for L = 4, 2, 0 are shown in Fig. 9. It is clearly seen that the two individual boojums coalesce into one large boojum as the distance between two vortex strings becomes small.
Let us next consider the asymptotically flat domain wall on which a vortex strings end from both sides. The moduli matrix is (4.5) The parameters z = ±L correspond to positions of the string endpoints. We show a typical solution with L = 6 form = 1 case in Fig. 10. The six panels in Fig. 10 show the same quantities plotted in the panels in Fig. 8. The magnetic force lines incoming from the vortex string on the positive x 3 side go into the other vortex string on the negative x 3 side. Seen from the x 3 axis, the distribution of the magnetic force lines is like a magnetic dipole, see the panel (a2) of Fig. 10. This is a remarkable contrast to the two vortex strings ending on the domain wall from one side. The former can be thought of as the dipole with opposite magnetic charges and the latter as the dipole with the same magnetic charges. As we reduce the separation L of the two strings, the region in which the magnetic flux expands gets smaller, see Fig. 11 where we plot the numerical solutions for L = 4, 2, 0. Note that when L = 0 the two vortex strings become collinear and the magnetic charges in the 2 + 1 dimensional sense exactly cancel. Even in the limit L → 0, the boojums do not disappear because the domain wall has the finite width and the boojums are separate in x 3 direction. Figure 10: The plots show the energy density isosurfaces of two vortices ending on one wall from two sides (a1, a2), where the blue and the red curves show magnetic fluxes, the wall energy density (b1), the vortex energy density (b2), the boojum energy density (b3) and the total energy density (b4) with the distance between two vortices L = 6.  Let us make a comment on the shapes of the boojums. When they are well separated, the shape of the individual boojum is approximately the same as that drawn in the right panel of Fig. 4. When the vortex strings are closer to each other than the vortex string size, the boojums merge. One may recall that rich 3D structure appears when several BPS 't Hooft-Polyakov monopoles in SU (2) gauge theory get close. For the boojums in the Abelian-Higgs theory, such drastic change in the shape is not observed, see Figs. 12 and 13.
Our last comment here is about the question of the definition of the binding energy raised in [65]. Consider two vortex-strings attached to the domain wall from both sides, see Fig. 14. Two vortex-strings are separated by the distance d. For this configuration, the authors of Ref. [65] argued that there are two possibilities where the binding energy is located. The first possibility is that the binding energy only localizes around the junction points. The second one is that it localizes not only around the junction points but also near the origin (dotted-circled domain in Figure 14: Schematic picture of the configuration of two vortices attached to the wall from both sided. the right figure in Fig. 14). Our numerical solutions, for example Fig. 10, strongly supports that the former is true.
Lastly, in Appendix A we present more complicated configurations with two or more domain walls.

Approximate solutions to 1/4 BPS Abelian master equations
In this paper, we have analyzed various 1/4 BPS configurations by numerically solving the corresponding master equations. However, in doing so, the need for asymptotically well-behaving initial configurations led us to the novel method, how to approximate composite solitons. In this section, we will collect our findings and present a general formula for an arbitrary configuration of vortex-strings interacting with up to two domain walls. The extension of this formula to three or more domain walls is, however, a straightforward task. Let us illustrate our method by recalling the simplest 1/4 BPS configuration, that is a vortex string attached to the domain wall, which we have studied in Sec. 3. In moduli matrix formalism, the corresponding master equation in dimensionless coordinates reads In Sec. 3 we argued that an approximate solution, which asymptotically approaches the true solution in all directions is is a solution to a domain wall part and u S (ρ) a solution to the vortex string part of the composite soliton. In other words, these functions are solutions to 1/2 BPS master equations The function U is designed to solve Eq. (5.1) in regions far away from vortex string, where |∂ ρ U| 1. Indeed, using asymptotic behavior of the string solution u S (ρ) ∼ log ρ 2 we can easily see that the function This shows that region of discrepancy between the true solution and U ∞ is not bound to a finite region like in the case of U ∞ . Amongst other advantages, the approximation U also enabled us in Sec. 3 to calculate all topological charges exactly and helped us to refine the arguments, which lead to the derivation of the general formula for the boojum charge in Eq. (3.51).
More importantly, the same strategy can be used to create approximations to virtually any 1/4 BPS configuration. For example, a collinear vortex strings attached from both sides of the domain wall is approximated as: Again, the reasoning is the same. Far away from the strings the ∂ ρ terms can be neglected and we are left with the master equation for a single domain wall, which is solved as On the other hand, far away from the domain wall |x 3 | 1 we can neglect ∂ 3 derivatives and solve the master equation approximately as u ≈ log em x 3 + e −mx 3 + u S (ρ) (5.11) with the first term on the right-hand side being an asymptotic form of the domain wall solution u W (x 3 ). The global approximation for two collinear vortex strings attaching the domain wall from both sides is thus u ≈ U = u W (x 3 ) + u S (ρ) as claimed.
Let us now provide a short list of other 1/4 BPS configurations and their global approximations. Their derivation follows the same arguments as presented above and we will not repeat them for brevity.
A semi-local vortex string of the size a attached to the wall is approximated as: where u SLS (a) is a solution to semi-local vortex master equation Two semi-local vortex strings of sizes a 1 and a 2 attached from both sides can be approximated as A general formula for approximate solution for arbitrary configuration of n 1 vortices attached from the right and n 2 vortices attached from the left is where the string solutions obey Generalizing further to a configuration of two domain walls with n 1 vortex strings attached to the right wall, n 2 vortex strings in the middle and n 3 vortex strings attached to the left wall, we find the approximation in the form: Here, we must supply double wall solution u W (x 3 ; δ) Arbitrary configuration of semi-local vortex strings attached from the right and arbitrary configuration of vortex strings attached from the left is given as And lastly, vortex string attached to a domain wall under the angle tan α = 2η/m is approximated as This ends our short list of approximate solutions, which we used at various places in this paper and in the second paper [67]. As we see, the general way how to construct approximations to composite solitons is to first solve the corresponding domain wall part u W by ignoring x 1 and x 2 derivatives and then replace these coordinates inside u W by appropriate vortex string solutions. This method can be thus applied to even more complicated configurations not mentioned in this section. As we see, to obtain fully analytical approximations we need to supply 1/2 BPS ingredients into our global approximations. We develop such analytic approximations for a single string and a single domain wall in the Appendix B.

Exact solutions to 1/4 BPS solitons
In this section, we will complement both numerical and approximative analysis by several exact solutions of 1/2 and 1/4 BPS solitons which we have found. Any composite soliton can be reduced to pure wall(s) or pure vortex-string(s) in some appropriate limit, i.e. either by shifting wall or string to infinity via the corresponding moduli parameter(s). Therefore, any exact 1/4 BPS solution potentially contains 1/2 BPS solutions as its limiting cases. Hence, before we can discuss exact composite soliton solutions we must first discuss exact domain wall and vortex string solutions.
In the following subsections, we will present a novel exact solution to a particular configuration of semi-local vortex-strings and an exact solution of a particular configuration of domain walls, out of which we construct new 1/4 BPS solutions. Throughout this section, we use dimensionful coordinates.
The master equation for the model with N F = 4 and the mass matrix given as above reads in general h 2 , h 3 , h 4 ) are elements of the moduli matrix. The generic configuration consists of a pair of domain walls, where the right wall has vortices attached from positive direction of x 3 axis (complex zeros of h 1 ), left wall has semi-local vortices attached from negative direction of x 3 axis (complex zeros of h 3 and/or h 4 ) and both walls have vortices stretched between them (complex zeros of h 2 ). The walls are bent if the net vorticity of incoming and outgoing vortex-strings is not equal. Asymptotically, the right wall has the profile x 3 ≈ 1 m log |h 1 | / |h 2 | , while the left wall has x 3 ≈ − 1 2m log |h 2 | 2 /(|h 3 | 2 + |h 4 | 2 ) . Let us first discuss an exact vortex string solution of this model. Since there are no known exact ANO vortices we set h 1 = h 2 = 0 and look for exact solutions over the leftmost degenerate vacuum. Getting rid of the x 3 dependence by shifting u → u − 2mx 3 in the master equation (6.1), we are left with the master equation for semi-local vortices: 2 This equation has an exact solution, namely for which |h 3 | 2 = 8 g 2 v 2 |z| 2 and |h 4 | 2 = |z| 4 . Loosely, we could think about u 2F SLS as a (non-linear) superposition of a local vortex and a semi-local vortex of the core size √ 8/(gv). The vorticity of u 2F SLS is indeed 2. Now, let us consider domain wall solutions. We set |h 1 | 2 = 1, |h 2 | 2 = e 2R and |h 3 | 2 = |h 4 | 2 = 1. The master equation (6.1) reduces to where the parameter 2R/m can be interpreted as the separation of walls for R m. If R is close to or less than m, the domain walls cannot be distinguished from one another and they form a single (compressed) domain wall. In the limit R → −∞ the middle vacuum disappears completely and only a single elementary wall remains.
There are two exact solutions for this equation. The first one was originally reported in [66]. Notice that u DW describes arbitrarily separated/compressed pair of walls (depending on the value of parameter R), but the ratio m/(gv) is fixed. On the other hand, the second solution describes walls with separations up to R = 1 2 log 2 − 4m 2 g 2 v 2 ≤ log √ 2. In other words, walls in u CW are always compressed, but unlike in the previous solution, ratio m/(gv) can be varied up to the upper bound 1/ √ 2. Reaching this limit, we have R → −∞ and the solution u CW reduces to an exact single wall solution, reported in [66], which we have denoted as u 1 in Appendix B.
Now we have all the necessary pieces to discuss an exact composite soliton solution, which reads with The interpretation of this solution is most explicit if we rewrite it as Remarkably, this formulation is functionally identical to the analytic approximation (5.2) for a single vortex string attached to the domain wall (only this time, the string is attached from the left side). 4 The bending of the compressed domain wall can be exactly obtained by equating both terms inside the logarithm in the solution (6.7), that is x 3 = − 2 gv log |z| 2 + 4 g 2 v 2 . Curiously, this bending is equivalent to the bending produced by a single semi-local vortex of the core size 4 g 2 v 2 . However, by comparing the solution (6.7) with the Eq. (6.3) we recognize that the attached string is identical to u 2F SLS . (The solutions become identical as x 3 → −∞.) Thus, the solution u 4F SW indeed describes a non-elementary wall u CW with nonelementary string u 2F SLS attached from the left. The same configuration with a string attached from the right side can be achieved by reflection x 3 → −x 3 . Notice that the position of the domain wall can be changed by shifting x 3 , while the position of the string can be changed by shifting z → z − z 0 . We plot the total energy density and boojum charge density on Fig. 15.
The master equation for the model with N F = 6 case and with the mass matrix given as above reads in general where H 0 = (h 1 , h 2 , h 3 , h 4 , h 5 , h 6 ) are elements of the moduli matrix. The generic soliton configuration is composed of four domain walls. As in the previous example, the complex zeros of the moduli elements h i (z), i = 1, . . . , 6 determine number and positions of vortices, attached and stretched between successive walls (as well as their asymptotic profiles).
Since we have doubly degenerate vacuum between the second and third wall, the exact vortex string solution of this model is again u 2F SLS , as defined in Eq. (6.3). However, it is necessary to introduce a more general solution of the semi-local vortex string based on N F = 3 model, with the master equation (6.11) That solution is u 3F SLS (s) = 2 log |z| 2 + s 2 , (6.12) with The solution u 3F SLS is well-defined only if the core size s is sufficiently big, namely s ≥ 2 gv . Indeed, if s = 2 gv this solution reduces to the previous case u 3F SLS 2 gv = u 2F SLS and bellow this limit |h 3 | is imaginary. The interpretation of u 3F SLS (s) is best seen if we rewrite the master equation as from which we see that u 3F SLS describes a pair of coincident semi-local vortices with core sizes s 2 ± 2s gv .
Let us now study the domain wall solution for the master equation (6.10), which is given as Again, the parameter 2R/m can be interpreted as the separation of walls for R m. However, we can rewrite u DCW as where e R = 2 cosh(mS). (6.19) In other words, the solution u DCW describes a pair of compressed walls u CW with separation 2S. 5 Now we have all the pieces to study an exact 1/4 BPS solution of Eq. (6.10). This solution is a combination of the wall solution u DCW with the semi-local string solution u 3F SLS stretched between its two domain walls: u 6F SW = 2 log e mx 3 + a + b |z| 2 + e −mx 3 , with the moduli elements given as where we have denoted The profile of the right compressed wall is asymptotically equal to x 3 = 1 m log a + b |z| 2 , while that of the left wall is x 3 = − 1 m log a + b |z| 2 . Therefore, the minimal 5 Notice that the solution is defined for all values of R, but the parameter S is only defined for R ≥ log 2. If R < log 2, it is possible to ascribe S a purely imaginary value S = i m arccos(e R /2), with which the solution (6.18) is still real. This confirms the intuition that the pair of wall merges at R = log 2, where separation becomes zero S = 0. separation between both composite walls is roughly 2 m log a. As in the previous case, we can rewrite this solution using its 1/2 BPS constituents: This is functionally identical to the approximate solution for a similar configuration given in Eq. (5.20).
Notice that the ratio of mass and gauge coupling (parameter α) can only vary in a certain interval, outside of which the solution is not valid. The edges of the interval correspond to extreme cases. In the limit α 2 → 1/2 we see that a, b → ∞, a/b → 4 and the solution becomes u 6F SW → u 2F SLS . In other words, this limit corresponds to the infinite separation of walls. In the other extreme α 2 → 1 we have a → √ 6, b → 0 and h 2 , h 3 , h 4 , h 5 → 0. As a consequence, both domain walls are infinitely compressed and the vortex string disappears. This domain wall is an exact solution, which we denoted as u 2 is Sec. 5.
By shifting x 3 and z coordinates we can change the position of the configuration as a whole, but interestingly the separation between the walls is not a modulus of the solution, but rather it is controlled by parameters of the model. As we can see, as the walls are further apart the core size of the vortex string 1 gv 2 1−α 2 shrinks ultimately to 2 gv . But, when the wall approaches each other, the size of the semi-local string fattens and, ultimately, diverges. This is a direct confirmation of a numerical observation which we make in Appendix A.
We plot the total energy density and boojum charge density on Fig. 16.

Other exact solutions
The two exact 1/4 BPS solitons u 4F SW and u 6F SW are just illustrative examples of the rich hierarchy of exact solutions. We believe that there exists 1/4 BPS exact solutions in multi-flavour models. Our preliminary investigations show that many new configurations of domain walls and vortices can be constructed as a suitable combination of its 1/2 BPS constituents, in the same spirit as described in this section. The list is potentially infinite if the number of flavors can be extended indefinitely.
To provide an exhaustive list of exact 1/4 BPS solution up to a certain number of flavors seems to be an important task, which however lies somewhat outside the scope of this paper. Therefore, we plan to address this issue fully in a separate work. For parallel domain walls, this was already achieved by one of us [69].
All our exact solutions presented in this section shares the same trait: they are only a special configuration of solitons. The domain wall solution u CW describes a pair of compressed elementary domain walls and both vortex string solutions u 2,3F SLS describes a pair of coincident semi-local vortices. In other words, some moduli are missing and we cannot provide exact results for generic configuration of solitons. Composite solitons u 4F SW and u 6F SW naturally inherit this trait. As a closing comment of this section, the existence of exact solutions of Abelian master equations directly implies the existence of exact solutions in non-Abelian case. This can be seen from the fact, that non-Abelian analog of the master equation can be, for special U (1)-factorizable solutions, decomposed into a set of independent Abelian master equations [31]. Thus, any exact solution in Abelian theory can be used to construct exact solutions in non-Abelian theories.

Outlook
In this paper, we have studied 1/4 BPS equations in the Abelian-Higgs theory. As written at the end of Sec. 6, our solutions also solve non-Abelian BPS equations by trivial embedding. While these embedding solutions are essentially Abelian solutions, genuine non-Abelian objects, like 't Hooft-Polyakov monopoles, may join the game in the non-Abelian theories. We will investigate non-Abelian extensions of 1/4 BPS solutions including monopoles, boojums, domain walls and vortex strings in forthcoming works. It is also interesting to study the 1/4 BPS equations in different spacetime dimensions. For example, the 1/4 BPS equations for vortex sheets and instanton particles in 5 dimensions [21,51], and those for the domain wall junctions in 3 dimensions [42,43] are known, but only a little is known about their solutions both in the Abelian and non-Abelian theories for the finite gauge coupling constants. We will also proceed to obtain numerical/analytical solutions for these series of 1/4 BPS equations. Furthermore, the technique developed in this paper may help to solve another type of 1/4 BPS equations or 1/8 BPS equations [70,53] for which nothing has been known except for the BPS equations.
We have stumbled upon unexpected analytic results, namely exact solutions to 1/2 and 1/4 BPS master equations in models with N F ≤ 6 described in Sec. 6. As we stated, we believe that these solutions are just simplest examples of a potentially inexhaustible wealth of exact solutions in models with an ever larger number of flavors of fundamental Higgs fields. The indirect evidence for this stems from the fact that our exact composite solitons, denoted as u 4F SW and u 6F SW , have been 'build' from of 1/2 BPS exact solutions for walls and strings. Therefore, it is natural to expect that more complicated combinations should give us another 1/4 BPS solutions in higher flavor models. The exact rules of this 'solitonic engineering' are so far unclear. However, the preliminary results show that certain restrictions on 1/2 BPS components must be in place if they are to be used in the construction of 1/4 BPS solutions. We have spotted glimpses of this in Sec. 6, namely that domain walls must have sufficiently low tension and cores of semi-local strings must be sufficiently large. The full analysis of this is a direction of future research. Furthermore, for parallel 1/2 BPS domain walls, such analysis was done by one of us in the study [69], where it was shown that exact multi-wall solutions can be build up from certain single wall solutions in quite an arbitrary fashion, except that the total tension of the final domain wall configuration must be small enough. Effectively, this restriction prevents the development of an inner layer inside domain walls, where unbroken phase appears. In fact, no exact solution of domain wall with this inner layer is known at present. Finding such a solution remains as an interesting open problem. The similar study for semi-local vortices is in preparation. By combining the findings for exact 1/2 BPS solitons we may be able to clarify, how their combination can yield exact 1/4 BPS solitons. Naturally, in doing so we may discover similar rules for other composite solitons, such as wall-wall and string-string junctions, which we did not address in this paper. An interesting question is, whether all exact 1/4 BPS solution can be decomposed as combinations of exact 1/2 BPS solutions, or whether some 'irreducible' solutions exist.
A. Non-axially symmetric solutions for N F ≥ 3 In this appendix, we describe configurations of two or more domain walls with a various number of vortex-strings attached. Throughout this section, we use dimensionless units.

A.1 Non-degenerate masses
Let us explain how to obtain the 1/4 BPS configurations with two domain walls separating three different vacua A (A = 1, 2, 3) in which n A vortex strings exist. The minimal model is N F = 3 with non-degenerate mass matrixM = (m 1 /2, (m 2 − m 1 )/2, −m 2 /2). The tension of the domain wall interpolating 1 and 2 is v 2 (2m 1 − m 2 )/2 and that of the domain wall interpolating 2 and 3 is v 2 (2m 2 −m 1 )/2. For simplicity, we will focus on the casem 1 =m 2 =m, namely that the two domain walls have the same tension.
First of all, let us consider the configuration with no vortex strings. The moduli matrix is characterized by only one real parameter δ as When δ 1 the physical meaning of this parameter is related to the separation of walls, since the position of the walls can be estimated as Thus the separation R is see Fig. 17. When δ is close to or smaller than 1, two domain walls coalesce and 2(log δ 2 )/m can no longer be understood as the distance. When δ → −∞, the two walls collapse into one heavier domain wall. The solution is parametrized as u = u W (x 3 ; δ) which satisfies the following reduced master equation Next, we put n A vortex strings in the A vacuum. The moduli matrix for this is given by H 0 = (P n 1 (z), δP n 2 (z), P n 3 (z)), (A.5) where P n A (z) stands for a monomial of n A -th degree. The corresponding gradient flow equation is One can rewrite this as Comparing this with Eq. (A.4), one finds a suitable initial configuration for solving the three-dimensional master equation Validity of this initial configuration is ensured by the asymptotic behavior u (n) S → log |P n | 2 for ρ 1. The positions of the two domain walls are estimated as Then the separation of the two domain walls is The domain walls are asymptotically parallel when 2n 2 = n 1 + n 3 and are asymptotically flat only when n 1 = n 2 = n 3 . In Fig. 18, we show the numerical solution for n 1 = n 3 = 0 and n 2 = 1 with δ = e 5 for the model withm = 2. Namely, we put a single vortex string in the middle vacuum 2 . Reasons for takingm = 2 here is, first, that the mass of each domain wall ism/2 = 1, and, second, that we have the exact solution u W to the domain wall master equation in this case [66]. Having the exact solution benefits our numerical works. We will mention more about this in Sec. 5. The panel (b2) of Fig. 18 clearly shows that the vortex string of the finite length = 2 m log δ 2 10 exists in the middle vacuum. The string length is the same as the domain wall distance with no vortex strings. The magnetic fluxes expressed by the blue curve in the panels (a1) and (a2) in the domain wall at the negative side of x 3 are squeezed and run through the vortex string, and then expand inside the domain wall at the positive side of x 3 . The endpoints of the vortex strings on the domain walls are accompanied with the boojums as shown in the panel (b3).
As δ gets small, the domain walls get close. At the same time, the distance between the boojums becomes smaller and the vortex string gets shorter. We show the solutions for δ = e 4,2,0 in Fig. 19. When the domain walls start to merge it is difficult to find a clear boundary between them. In this situation, it is difficult to speak about the length of the string. However, the diameter of the string becomes very wide, while the magnitude of the magnetic flux weakens, see Fig. 20 for δ = e −2 . The panels (b2) and (b3) in Fig. 20 show the vortex and boojum charges and we see that the vortex is indeed sandwiched by the domain walls. Ultimately, at δ = 0, the middle vacuum completely disappears, so that the vortex vanishes as well. Fig. 21 shows changes in the shape of two boojums living in the different domain walls. As the domain walls close on each other, the boojums merge into a flat disk. Figure 18: The plots show the energy density isosurfaces of one vortex stretching between two walls (a1, a2), where the blue and the red curves show magnetic fluxes, the wall energy density (b1), the vortex energy density (b2), the boojum energy density (b3) and the total energy density (b4) with the half distance between two vortices R = 5. Figure 19: The plots show the energy density isosurfaces of one vortex ending on two walls. The half distance between two vortices is taken to be R = 4, 2, 0 from top to bottom. Figure 20: The plots show the energy density isosurfaces of two vortices ending on one wall from two sides (a1), where the blue and the red curves show magnetic fluxes, the vortex energy density (b2), the boojum energy density (b3) and the total energy density (b4) with δ = e −2 . Figure 21: The plots of the boojum energy density for the case that two vortices ending on one wall from one side (top), two vortices ending on one wall from two sides (middle) and one vortex ending on two walls (bottom). The half-distance of two boojums is R = 4, 2, 0 from left to right. As an example for more complicated configuration, in Fig. 22, we show the solution which have single, double and triple vortex strings in the vacuum 1 , 2 and 3 , respectively. The number of vortex strings satisfies the relation 2n 2 = n 1 +n 3 , so that the domain walls are logarithmically bending but are asymptotically parallel.

A.2 Partially degenerate masses
When some masses are partially degenerate, the corresponding vacua are degenerate. Then the vortex strings there become semi-local vortex strings, which we discuss in the second paper [67]. As before, we consider the minimal model N F = 3 with M = diag(m/2,m/2, −m/2). The relevant moduli matrix is given by where P n (z) is a monomial of power n and Q n (z) stands for a polynomial of power n. The n 1 semi-local vortex strings in the degenerate vacuum 1 are determined by P n 1 and Q n 1 −1 , while the n 2 local vortex strings in the vacuum 2 are determined by P n 2 . The gradient flow equation to be solved is This is identical to Eq. (4.2) if we replace |P n 1 | 2 → |P n 1 | 2 + |Q n 1 −1 | 2 . Thus, a suitable initial configuration is given by where u  SLS . These pieces can be supplied as numerical solutions of corresponding 1/2 BPS master equations, which are generally much easier to solve than the full three-dimensional 1/4 BPS master equations. A complementary way is to obtain analytic approximations to 1/2 BPS solitons. This is the goal of this Appendix. In the following subsections, we will develop accurate analytic approximations to both vortex string and domain wall.
Especially accurate approximations can be found for vortex string master equation Eq. (5.4), since there is no intrinsic mass scale as in the case of domain wall master equation Eq. (5.3). This allows easy use of the so-called global Padé method, which we describe in the next subsection.

B.1 Approximations to ANO vortex
The master equation for single vortex string reads No analytic solution of Eq. (B.1) is known. In order to construct global approximations, we take the advantage from what we can learn about u S locally. In particular, let us investigate behaviour of u S close to ρ = 0 and for ρ 1. Expanding u S in Taylor series around origin and matching coefficients on both sides of Eq. (B.1) we obtain where u 0 ≡ u S (0). Note that have used the regularity condition u S (0) = 0. On the other hand, the asymptotic behaviour of u S is u S (ρ) ∼ 2 log(ρ) + qK 0 (ρ) ∼ log ρ 2 + qρ 3/2 e −ρ as ρ → ∞ .
Here, K 0 (ρ) is the modified Bessel function and q > 0 is some constant. Both u 0 and q are locally undetermined; their values cannot be ascertained neither by Taylor series nor by asymptotic series. In other words, these numbers characterize global properties of u S and as such, they are very hard to study analytically. Recently, a detailed study of Eq. (B.1) was published [71], where based on a perturbative expansion around a small winding number, the numbers u 0 and q (2D 1 and 2C 1 in their notation) were obtained with high precision, both analytically and numerically. In the following, we shall adopt the values u 0 = 1.01072165 and q = 3.41572835 (Eq. (2.66) in [71]) in all our approximations.
A (global) Padé approximant is a widely used tool to model a function based on its series expansion(s). The basic idea is to find a rational function which, when expanded, match all series to given orders. Generally, such approximant tends to be very accurate even outside radii of convergence of individual series and for sufficiently nice functions it remains close to the true solution everywhere.
Nevertheless, there is only so much a rational function can do. In particular, as we see from Eq. (B.3), the asymptotic formula for u S contains log, exp and square root (suggesting essential singularity at infinity), neither of which can be globally approximated by rational functions. Therefore, simple Padé approximants are insufficient in this case.
A way to overcome these difficulties, originally discussed in [72], is to put the Padé approximant directly into the asymptotic series, in such a way that singular terms are replaced by rational functions, which coefficients are fixed to match Taylor series at the origin. In this way, the resulting function is guaranteed to behave in a desired way both at the origin and at infinity.
Let us, therefore, propose the following ansatz: where P p+3 (ρ) and Q p (ρ) are polynomials of order p + 3 and p respectively. The coefficients are chosen such that they give correct large ρ behaviour and the remaining ones are fixed to match the Taylor series (B.2) up to O(ρ 2p+3 ). After this procedure is done, it needs to be checked that both P p+3 (ρ) and Q p (ρ) are strictly positive for ρ ≥ 0, as we require that u (p+3,p) is regular on the positive semiaxis. It turns out that this is satisfied only up to p = 4 and for higher p singularities appear.
We have found the following Padé approximants  Let us illustrate the accuracy of our approximations for p = 0, 2, 4. In Fig. 25 we show relative difference between u (p+3,p) and numerical solution and in Fig. 26 we show the master equation evaluated for each approximation. It demonstrates that u (p+3,p) is increasingly more accurate for larger p. As we already mentioned u (8,5) is not a regular function. This seems to persist even for higher values of p, indicating that capacity of this particular ansatz to approximate true solution has been exhausted.
Approximations u (p+3,p) have, however, one disadvantage, namely that they depend on odd powers of ρ. In contrast, series (B.2) contains only even powers. The coefficients must be therefore fine-tuned to eliminate all odd powers, especially the first power, since its presence cause a singularity in the master equation (B.1) at ρ = 0. Due to the finite precision of decimal representation of real numbers, this fine tuning is realized only imperfectly, which means that all odd powers, although strongly suppressed, are still present. (dashed yellow line), u (7,4) (red dotted line) and uc (6,3) (purple dash-dotted line). Notice that we had to scale down the first graph by a factor 9000 and the second graph by 20 to fit them comfortably into the picture.

B.2 Approximations to domain wall
In this subsection we shall for notational convenience relabel the third coordinate as y ≡ x 3 and we also omit the˜on the parameter m ≡m, even though we are still using the dimensionless coordinates. which gives δu W = c e −y + 1 1−4m 2 e −2my , where c is an arbitrary constant (the second independent solution to homogenous equation e y is discarded, since δu W must be negligible compared with my for large values of y). This formula does not apply in the case m = 1/2, where instead we have δu W (m = 1/2) = c e −y + 1 4 e −y (1 + 2y). As this case is somewhat anomalous, we shall omit it and in the following, we assume that m > 1/2.
Combining the δu W correction with the dominant term into a single logarithm we obtain the asymptotic form At this point, we do not know the value of the constant c. However, for specific values of m, one can exploit reflection symmetry of the solution u W (−y) = u W (y) to obtain a rather good estimate. In order to do that, let us look on additional sub-leading terms. Repeating the procedure outlined above, we obtain next order corrections combined into a single logarithm as: Notice that this asymptotic form contains (upon factoring out the dominant term e my ) either powers of e −2my or powers of e −y or their combinations. If we focus on half-integer values of m these factors can be the same, which we use to our advantage, since in that case, we can force the reflection symmetry. Let us illustrate this on the case m = 1. We have u 1 ∼ log e y + c + c 2 − 1 3 e −y − c 6 e −2y + 4 135 e −3y . (B.29) The strategy of harvesting approximate solutions from the asymptotic forms, which we adopt from now on, is that we discard terms not compatible with the boundary condition u → −my as y → −∞. In the present case, the terms e −2y and e −3y would grow faster than e −y as y → −∞, so we discard them. Doing so it only remains to ensure the reflection symmetry of u 1 , which fix the constant to be c = 2. Thus u 1 ∼ log e y + 2 + e −y . (B.30) Surprisingly, this approximation is, in fact, an exact solution (B.23)! It turns out that terms we have discarded are exactly canceled by sub-leading terms. The next case m = 3/2, including next-to-next leading order corrections to u, has the asymptotic form with incompatible terms already discarded given as and unique value of c, which ensures reflection symmetry, is now c = 2 √ 6. Indeed, this choice of c makes u 2 an exact solution (B.25).
All three exact solutions u 1 , u 3/2 and u 2 were first reported in [66] and it seems that they are the only ones. Indeed, if we continue this procedure for m = 5/2 we obtain the asymptotic form u 5/2 ∼ log e 5y/2 + c e 3y/2 + c 2 3 e y/2 + c 3 24 e −y/2 + c 4 540 e −3y/2 + c 5 − 1080 25920 e −5y/2 (B. 33) and upon inspection, we conclude there is no choice of the constant c, which would make the above approximation an even function of y. At this point, we adopt the attitude of 'patching' the asymptotic form in such a way, that it becomes symmetric under the reflection y → −y. The condition for c is chosen in such a way that the patching affects only the innermost terms. Thus, we choose c 2 /3 = c 3 /24, which yields c = 8. Then we fix the coefficients of e −3y/2 and e −5y/2 to match their reflected counterparts. In other words, we have u 5/2 ∼ log e 5y/2 + 8e 3y/2 + 64 3 e y/2 + 64 3 e −y/2 + 8e −3y/2 + e −5/2 , (B.34) Figure 27: Relative errors between numerical and approximate solutions u 5/2 , u 3 , u 7/2 and